library(readr)
package 㤼㸱readr㤼㸲 was built under R version 4.0.5replacing previous import 㤼㸱lifecycle::last_warnings㤼㸲 by 㤼㸱rlang::last_warnings㤼㸲 when loading 㤼㸱hms㤼㸲replacing previous import 㤼㸱lifecycle::last_warnings㤼㸲 by 㤼㸱rlang::last_warnings㤼㸲 when loading 㤼㸱tibble㤼㸲replacing previous import 㤼㸱lifecycle::last_warnings㤼㸲 by 㤼㸱rlang::last_warnings㤼㸲 when loading 㤼㸱pillar㤼㸲
library(ggplot2)
library(dplyr)
package 㤼㸱dplyr㤼㸲 was built under R version 4.0.5
library(rms)
package 㤼㸱rms㤼㸲 was built under R version 4.0.5package 㤼㸱Hmisc㤼㸲 was built under R version 4.0.5package 㤼㸱SparseM㤼㸲 was built under R version 4.0.4
library(synthpop)
package 㤼㸱synthpop㤼㸲 was built under R version 4.0.5
library(tidyr)
package 㤼㸱tidyr㤼㸲 was built under R version 4.0.5
library(plotly)
package 㤼㸱plotly㤼㸲 was built under R version 4.0.5
library(sjPlot)
package 㤼㸱sjPlot㤼㸲 was built under R version 4.0.5
library(performance)
package 㤼㸱performance㤼㸲 was built under R version 4.0.5
# install the 'mi' package if you are planning on redoing imputation on the original dataset. For now the imputed dataset has already been saved
# library(mi)
source('./helper functions.R')

Reading in original dataset, performing multiple imputation and writing out imputed dataset. This chunk commented out because we’ve already done this and saved the resulting file as our starting point.

Read in dataset from Brinati paper. Dataset already been imputed based on code above.


brinati_covid_study_v2.imputed <- read_csv("data/brinati-covid_study_v2_imputed.csv",
  col_types = cols(GENDER = col_factor(levels = c("M",
  "F")), SWAB = col_factor(levels = c("0",
  "1"))))

glm.orig.fit<-glm(SWAB ~ ., brinati_covid_study_v2.imputed, family='binomial')
glm.fit: fitted probabilities numerically 0 or 1 occurred
glm.orig.fit

Call:  glm(formula = SWAB ~ ., family = "binomial", data = brinati_covid_study_v2.imputed)

Coefficients:
(Intercept)      GENDERF          AGE          WBC    Platelets  Neutrophils  Lymphocytes    Monocytes  Eosinophils    Basophils          CRP          AST          ALT          ALP          GGT          LDH  
 -2.185e+00   -8.505e-01    1.746e-02   -1.146e-01    2.011e-03   -1.289e-01    9.374e-03    8.215e-01   -1.095e-08   -4.520e+00    5.427e-04   -7.082e-03    1.026e-02   -1.195e-02    5.287e-03    1.001e-02  

Degrees of Freedom: 278 Total (i.e. Null);  263 Residual
Null Deviance:      366.4 
Residual Deviance: 239.4    AIC: 271.4
write.csv(summary(glm.orig.fit)$coef, file='output/OR-original-model.csv', row.names = FALSE)
predicted = plogis(predict(glm.orig.fit))
observed = brinati_covid_study_v2.imputed$SWAB
assessPerf(predicted, observed)
          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p             U      U:Chi-sq           U:p             Q         Brier     Intercept         Slope          Emax           E90          Eavg           S:z           S:p 
 7.322477e-01  8.661239e-01  5.000567e-01  4.514185e-01  1.269458e+02            NA -7.168459e-03  8.242296e-13  1.000000e+00  4.585870e-01  1.361516e-01 -2.780966e-08  1.000000e+00  6.627578e-02  5.261796e-02  2.628540e-02 -2.750989e-01  7.832403e-01 
p<-plot_model(glm.orig.fit)
p<-ApplyFigureThemeLargeFontOnly(p + ylim(.1, 2.5) ) 
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
  
p
ggsave(filename = 'figs/OR-original-model.png', plot=p, width=6, height=6, units='in', dpi=300)

Draw ROC curve

cms.cutoffs <- lapply(seq(0.01,0.99,0.01), function(cutoff) {
  ret=confusionMatrix(predicted, observed, cutoff=cutoff)
  ret$cutoff = cutoff
  ret
  })
ff<-data.frame(t(sapply(cms.cutoffs, function(cf) c("Cutoff"=cf$cutoff, "Sensitivty"=cf$sens, "Specificity"=cf$spec))))
colnames(ff) <- c('Cutoff', 'Sensitivity','Specificity')


cutoffs_to_plot <- c(0.09, 0.3, 0.6, 0.9)
p <- ApplyFigureTheme(ff %>% mutate(fpr=1-Specificity) %>% filter(Cutoff %in% cutoffs_to_plot) %>%
  ggplot(., aes(x=fpr, y=Sensitivity)) + 
    geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
    xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
  )
ggsave(filename = 'figs/example-sens-spec-on-roc.png', plot=p, width=6, height=6, units='in', dpi=300)

for (cutoff in lapply(1:length(cutoffs_to_plot), function(i) cutoffs_to_plot[1:i])) {
  p <- ApplyFigureTheme(ff %>% mutate(fpr=1-Specificity) %>% filter(Cutoff %in% cutoff) %>%
    ggplot(., aes(x=fpr, y=Sensitivity)) + 
      geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
      xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
  )
  ggsave(filename = paste0('figs/example-sens-spec-on-roc-', cutoff[length(cutoff)],  '.png'), 
         plot=p, width=6, height=6, units='in', dpi=300)
}
p <- ApplyFigureTheme(
  ggplot(ff, aes(x=1-Specificity, y=Sensitivity)) + geom_line()+ 
    geom_ribbon(aes(ymin = 0, ymax=Sensitivity), color=NA, fill='blue',alpha=0.3) + 
    xlim(0,1) + ylim(0,1)+
      xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)') 
)
ggsave(filename = 'figs/example-roc.png', 
         plot=p, width=6, height=6, units='in', dpi=300)
p <- ff %>% mutate(fpr=1-Specificity) %>%
  ggplot(., aes(x=fpr, y=Sensitivity)) + geom_point(aes(frame=Cutoff)) +
    xlab('False Positive Rate (1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
  
  
ggplotly(p)

Draw calibration

```r
cal_breaks = seq(0, 1, .1)
cal_deciles <- data.frame(Predicted=predicted, Observed = as.integer(as.character(observed))) %>%
  mutate(Bins = cut(Predicted, breaks=cal_breaks, include.lowest = TRUE))
cal_deciles_summary <- cal_deciles %>% group_by(Bins) %>% summarise(n=n(), Average_Observed = mean(Observed),
                                                                    Average_Predicted=mean(Predicted)) 


for (cutoff in lapply(1:nrow(cal_deciles_summary), function(i) cal_deciles_summary$Bins[1:i])) {
  p<-ApplyFigureThemeCalCurvePoints(
    ggplot(cal_deciles_summary %>% filter(Bins %in% cutoff), aes(x=Average_Predicted, y=Average_Observed)) )
  p<-ApplyFigureThemeLargeFontOnly(p)
  SaveStdSquareFigure(p, filename = paste0('figs/example-cal-curve-points-', cutoff[length(cutoff)],  '.png'))
}

p<-ApplyFigureThemeCalCurvePoints(ggplot(cal_deciles_summary , aes(x=Average_Predicted, y=Average_Observed)) )
p<-ApplyFigureThemeLargeFontOnly(p)

SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all.png')

p<-p + geom_smooth(method='lm', se=FALSE)
SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all-with-bestfit.png')
p

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArwAAAGwCAMAAAB8TkaXAAAA8FBMVEUAAAAAADoAAE0AAGYAAP8AOjoAOmYAOpAAZpAAZrYzMzMzZv86AAA6OgA6Ojo6OmY6OpA6ZmY6ZpA6ZrY6kLY6kNtHR0dNTU1mAABmADpmOgBmOjpmZjpmZmZmZpBmkJBmkLZmkNtmtttmtv+QOgCQOjqQZgCQZjqQZmaQZpCQkGaQkLaQtpCQttuQtv+Q2/+2ZgC2Zjq2kDq2kGa2kJC2tpC2tra2ttu227a229u22/+2/9u2///bkDrbkGbbtmbbtpDbtrbb27bb29vb2//b/7bb/9vb///r6+v/tmb/25D/27b/29v//7b//9v////kkeULAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAdWElEQVR4nO2da2PbOHaGqbHHaWLNpI03ltdpumni1t3tpFE67bR1o062ma0jOZb+/78p7wKICwHiQhzyPR8SW5YfHoJPkCMQBLIDAkE0srETQCCGBuRFkA3IiyAbkBdBNiAvgmxAXgTZCCTvMwQiXASWV/Ozrx74PhiAkISsVivICwhFSK7u6gB5AaEHWZXuQl5A6EEqdSEvIOQgjbqQFxBikNXRXcgLCCkIoy7kBYQShFPXXt7Hq6f8Cw+3WZa9uJN8w+L7shoaabYxIEEgq4671vKuM17ex6usiMUH4RsOr89qeKTYxoCEgXTVtZV3v854efc32ZO7w0P+533nGx6vzcohEmxjQIJARHUt5X24zjry7pYnnw5ll/uu8w2P12TlFMm1MSBBIELFUIaNvNssu/iVl3dTf7vJXna+4fHKrBwjsTYGJAhErq6lvGfvc4E5edd1J7stSgXuGx6vyMo5kmpjQMJAFOraf2Dj5c2r3MrX3fLJPfdNBy/Pyj1SamNAgkAuVeq6yvt4VY8sFOUu982hnsv7FYFwiMsiVD8MKC/7b0MW6HkB6Yu8171U/xTyApIspCx2NRDUvICkATk/77xQjzEEkxejDYD4gZwXwb3SjDGEkxfjvID4gCjVDSkv7rAB4g6RqtuMj4WTN69zz9i5DWeY2wCIJURdMfRBBstbTy/bLZmJZNw3HF4WkBeQg7Zi6IO4yov5vIC4QLQVQx8ET1IAMiKkr9vVQyAvIKNBDNSFvIAkCem6K5/5CHkBSQ5i1O3qIZAXkFEgsm7XFgJ5ARkBYlYx9GUCeQGJDzHudvWZQF5AYkNs1IW8gKQE6birqRj6MoG8gESF2HW7+kwgLyARId1JOKy6p6enlplAXkDiQTQVw2kVVplAXkBiQbQVA+QFJF2IrmJo3ZXYC3kBGRvSM8YAeQFJFdI7PAZ5AUkUoq0YqkDNC0iKEAN1IS8gSUJ4d9U31DDOC0hiEKNud2AmkBeQkBCx2/WYCeQFJBzEtGIYmAnkBSQYxL3b1WcCeQEJBPGiLuQFZAQI5+6wiqEvE8gLSAhIqW4LGa4u5AUkMqSehFNDXNSFvIDEhTQVQwlxqBj6MoG8gHiGHIvdAuKoLuQFJB6Enbar2f/PSyaQFxCfEG6M4dKDu2PKG3PDOcTYUah7/E63/5+XQM8LiDdId2hXs/+fl0wgLyCeIOJdidCZQF5A/EAEd8NnAnkB8QGRqAt5e4Lu1Z4WhHX3eFcC8mqD7NWeFETa7UbIBPIC4ghhbkucnn5nuP+fl0wgLyBuEE7dIuJlAnkBcYGwFUPpLvsIMOTVBr2rPSkIt/5YpS77ADvk1Qa1qz0tCD/G0LgLeQ2D2NWeFKQ7PPadsNwY5NUGqas9LUjXXcmKTZBXG5Su9qQgorqQ1zboXO1pQRh3O2vzx8wE8gJiDZF1u2NkAnkBsYV0ut3xMoG8gNhBFBXDCJlAXkDsIMbdbvBMIC8gVhAbdSFvTyR/tacFObprtJQI5NVG6ld7UhC7bjdkJnVAXkDMIOd8tztiJm1A3qlDFPuU2GZiWTHIIUMC8s4Wotwhyi4T64pBBhkWkHe2EC/yDqgYRMjQgLxzhah3RbWADKkYBMjggLxzhXiQd1jFIGQyOCDvXCHu8g6sGIRMBgfknS3EseYdXDGImQwNyDtbiJu8Lt1uN5OhAXlnDBk+zuuoLuTtiWSVmQCkdXfwlihJyftwm2XZi7vjC/ubrImTT4fD49Xxax5vmZVxJHW1JwVx7na9ZaL+kY28tZuLD+0rHXl3S8g7EQjX7Y6aiR95c1Of3B0e8j/vuz/aLUujt9lTBd4yK+NI52pPCuJeMfjKxJe8u2XZo+b977vOT3KtXxZ/r6u/JHjLrIwjlas9LYiXbtdLJnqIhbybul/dCIpuqs54f8NUFDzeMivjSORqTwpSqFsOr4Xc/88LxELedd3jbrt1Q9MXP16d/HSdZWcfRLxlVsaRxNWeFOS8cfc7903U0pE3rw0qeXfLjrzr2ubm8xrTM0NeYpBW3TxC7v/nBWIu7+NVXRTUtS/zg7ZLLj7RfbutxyOeFRF0EzmE5yg3AWzcPR07m97wIO+mKSOaL9bHQQf0vIQg9bTdWt1BN+Y8ZWIIcZe3GWo4BlMUQ146kGaMYeXL3YTkVdW83Y6YewXyUoEww2PdBc4jZ2IOcR9t2Ah3JiAvPUjb7TYr7Y6WiQXEfZy3vTPRds0oG6hBWHXj7P/nBeJ8h60thdsPamwRDHkpQPhud8xM7CB2cxvOxLkNTI2wW2YX94eHa0zMIQWRqEvkdGxmldU3IaqethkPYyvgTdaddgZ5k4fU7nam4JA4ncHzeWXyHh5e5epeMB0z5E0cIu12R8lkAARPUswZcs50u+NmMggCeWcMkVcMY2QyDAJ5ZwtRVQzxMxkKgbwzhagrhtiZDIdA3nlCNBVD5EwcIJB3jhBtxRA1EycI5J0hRF8xxMzEDQJ5ZwfpqxjiZeIKgbxzgxh0u5EycYZA3nlBzNQlcjqQd1aQyl2DpURInA7knRHEtNsNn4kfCOSdD+TY7Y6diScI5J0LxLhiCJ6JNwjknQnEptsNm4k/COSdBcRS3dRPpw7IOwNINQnHarnSlE+nDcg7fYi9ukmfzjFEeb994eM3lyND3tEh1hVDsEy8QwR5m30luPX6BwfkHRlyrBjGziQABPJOGjKkYgiTSQhIr7zfP4e8VCHDKoYQmYSBKGref8wWb3877H9dSnZPsQnIOyJkYMUQIJNAEPlowzY7rhfd3YDCKiDvaJDBFYP3TIJBpPLub5pKNy8inLpeyDsWxKXb9ZtJOIhU3lzZo7z4wEYQ4qhuaqejCJW8dbWwwWgDRUjhbtr7/3mByGvedbEm2S+//Pur/O/u0tFWAXlHgDh3u94yCQyRy9tuSsUt+TggIG98SNPtjp/JSHMbfl16cRfy+oMYLrXvXjH0ZpIMRDUx59vPPxR3KN44jfJCXm+Q01OzTU68dLvaTBKChJ5VNvY+c1OJWt6edxXqXuYRJaXxQy3vtz//h2O3e0DP6wtyemrQ9Z5X3e6lc7eryyQliEreoug9+fT4+793OzLk9QMxkbepGAicjieIQt6P1Xyyxyu3kTLI6wnSL++x2CVwOp4gyrkNjbzcvlXWAXk9QXrcPWc+qFE4HT8Q5U2Kp3+5OvlU3GHD3IYUIHp5uTEGCqfjB6K6Pbz48FjIi7kNyUAM1aVyOj4g6ok5kJcIpDu0S/x0LCA9Pe8WZUPiEPGuBOnTsYIoa94n/3l18t/FgBkm5qQMqd1NIJMRIJiYQxkivRlM93RsIaqJOa28Tk8BQd6gEEm3O1Im40BUd9iKXYTzOHvvdmTIGw4iV5fs6QyAaCbmfPnyxfnIkDcU5FzlLs3TGQSRjzb88c7HYQ+QNxhEqS7N0xkGUQ2VOc/k5fGygLyDIRp1KZ7OUIjyAUz3epfFW2ZlHCTa2DekrhgSyGRkiHzdhn+ohxpeuJYPkNc/RNvtRs1kbIjiA9v+83U9zOtWPkBe35A+dYmdjhNEPdqwL59iw7oNaUH0FUPMTBKAaJ9h+3YLeZOC9He7sTJJAqLpeevKAfImAzHodiNlkgZEJe/nV6h5E4OYqUvmdDxA5KMNHzHakBzE1F0ip+MDgnFeGhBjdWmcjh+IcjL6G6ddgLp4WUBec0jpbhKZpARRyPs7H4c9QF4/EItuN3AmaUFUK6O7bUUh4GUBeQ0hNt1u2EwSg/SsjO4akNcZYqlu6qfjE6LqeR2XNu3iZQF5DSDn1u4mfTp+IYoVc5bZmZcpvZDXDWKvbtKn4xki73n/qbq5tvgxD2wiOBpkiLoJn453iGacF9u3jgsZUDEEyiRVCORNFTJQ3VRPJwREXjb8+ZdjMEtMP9wKt4wbzyvFhTdA3oGQweqmeTphIDbL+teicgMRzfIkpbziGyDvMIiDuymeTiCIxbL+xa2Lu8MDfwNjyywHJXkD5B0CcVE3wdMJBrFY1n+3bLpXZr3pNbOkjuQNkJePnl0lKki5CWDoTKYAsVjWf1N/vWF83d8wNYLkDZCXjd59JQqIW7drmMk0IBbL+q/rL7dMWfB4dfLTdZadfVC8AfKyYSKvY7drmMk0IObL+ucVbeXmbnmUt11O8qX8DZCXif5dUb66q0usTZwg5sv6P17VFUJd2pZRLD59Vzyomf+s+4ZnRYy9z1xK0cqresNl4W7MjIiH+bL+Unk3dce8zstd2RvQ8zLR0/OufHS7ZplMBGK+rL9U3iaKMhfy9kUUd4m1iQvEfFl/ac3bROErat6+0MhbqUvrdEaHWCzrLxttOP5C3tlitKE3dOqe0zudkSEWy/pLhnHbzrb0FeO8A6OtGEbPhBbEYll/2R22deVrLvFL3GEbGEyxO4XTiQixWNY/N/SsO7chry8u7g8P1+WYhOQNkLcvmoph/EzoQfSzyr5xxW1dCVdVcN3nbtjKmHsDh7fMyjhItLEu+DEG8qcTF6IsG/74qXiSjX+UjZ2uW8tb1heLi3vxDRzeMivjINHG6ugOjxE/ndgQhbybvA6o+lG3Z+AhryZWwtAu6dOJD1FOzDn5tPawiyDkVYfkrgTl0xkBorxJ8dd/ucqyJ1/+FRtnh4GI3e5YmdCFqG4Pvys/fb2T3gu2CMgrD6m6dE9nJIhC3mo2ZDlVDPL6hyjcpXo6Y0GU8u5vyjk56Hn9Q1TqEj2d8SCqtcreFlXDy8N+jZrXM2TF3ZYYMxPyEOUHtvJmQ/EYEEYbvELU3W7sTOhD5PJWCzA8dV/rFPLyoVWX3umMDFHcpHi4zopJC49XZ27r9EJeNnQVQ9xMpgHRz23Yf5G+bB6Ql4mebjdiJhOB2Cz3NCAgbxv96pI6nRQgSnkf/vT69WvnLYEgbx0rE3fpnE4aEE3NW20j6LaxCuStwkhdOqeTCEQ32lA9S4HRBmeIWbcbI5NpQTTjvM9fv8KsMg8QY3VpnE5CENUdtuyvinphv85wh80RYuEuhdNJCaJedKT8arfEsv5OEBt1CZxOWhBVz1srm2uMnnc4ZGXnbuqnkxpEXvNump53mzEL9A6IWctrq27ip5MeRC5v9RD74fBw5bgV5pzltXc36dNJECLI+/iq2DlwWW4hWPz5/O9QNgyBDFA35dNJEiLKy+7Bhn3YBkIadS2XfUz1dBKFQN4QkEHdbpBMJg0R5OU2EOxsImgfs5R3qLqJnk66EMwq8w0ZWDEEyGTyEMjrGTK82/WdyfQhqlll//b69et/vhPebhtzk9dJ3fROJ3GIVN7PP9Qf1r5/63jkecm7cnQ3sdNJHiKT9yMz1nDhduRZyeuqbmKnkz5EIu+GGylzmhE5J3kLdbFtcFSIKG+x2sjiTTE+9u3nZeZ6e3jsfeZixWUR2AMwbojybhhhq2VzXOTV/GxKPe9qdenc7frJZE4Q8SbFDb8DEKZE9kdRMWDb4PgQ2e1hplLAZPT+WPmodr1kMjeITF5GV9f1nmYgrz91kzgdShDI6wYp1fXl7vinQwsirXmPmwRuUfPqgqkYSFztiUHE0Qb2geFifiQeA1JGqy62DR4FIsqbd7bZk2pWQ7luzjvJbxnHlOXlKwYSV3tiEFHeom7IssXfvP7bZbVIr0tMV97uGAOJqz0xiOT28OM1c3cY6/PKg60Yxs1kxhDZxJz9beuu4zp7U5VXMsZA4mpPDCKfz/vwpx+X2eLHPzhP6J2kvNK7EiSu9sQgeJLCGiIf2iVxtScGgbyWENVdCRJXe2IQyGsFWSnvqJG42hODQF4LiFpdIld7YhDIaw6p1FVMZCBxtScGgbymEK26RK72xCCQ1wyy6nGXxtWeGATyGkH61CVytScGgbwGkFpd7bRdEld7YhDI2wvprRiiZQIIH5C3D8Koq3tagsTVnhgE8uohJhVDnEwAEQLy6iBmFUOMTACRBOTVQAwrhgiZACILyKuEGFcMwTMBRB6QVwGxqBgCZwKIKiCvHGLX7YbMBBBlQF4ZxEzd09PT4JkAAnmtICsjd0+rCJoJIJDXDmJYMUDe0SGQtwNp1DV0t7WXxNWeGMRK3ofimfgX/CPFn19l2aJ+rdk987g2HzV5zSqGIiDv+BAbeWs3uYX+681XFuWiULslcXlZdXsGGSDv+BALefc3xRpmDzfsupHbbPH2ULxW+roVFociJa9pxVAFat7RIRby7paloXn/2669l/tcbgFQv7YWNrAgJK95xVAF5B0dYiHvpu5XN0dFH6/qGqLUdn8jbB1ER16LiqEJjPOODLGQd133uFvJetOlvI9XJz9dZ9nZBxFvmZVx+Goeu4ohZCaAmEPM5c1LhEre3VKQt+qBm89rTPFAQ17biiFcJoDYQMzlbUuEuvZlY1N2xttyVepvt/V4xLMixt5nziSK/f/KL7AJIK3wIu+20nVT1xPr46ADgZ633P+vCGwbTA3iQ97tcsGt/c8UxcnLWxQMl+VX2DaYHMRDzbvp7k/M2J24vFWxW0CcN6MicbUnBnEebdh/FPbWJiNv/Tntq4+N1Ehc7YlB3MZ5y/642bai7ZqJlA3t8Bi2DSYKcbvDdtivmW64/qDW3HZj8ZZZGcdwRjs85mn/ShJXe2IQu7kNZ925DRv2m90yu7gvNm8jMDHHs7pErvbEIDazyuqbEFWJW3azzRzIer+2TdaddpaovEd1sW0wYcjg+bylvNuMk/fwUEzuvWA+ziUpb7diIHGhABFjhk9S8Opi22C6kNnJK6kYSFwoQMSYmbzSMQYSFwoQMSYlLzfBVhZCxeAtEUBGgExI3s6jDWKoxhhIXChAxJiPvOq7EiQuFCBiTEfe7uO8ndAM7ZK4UICIMRN5tXclSFwoQMSYhbwr/R01EhcKEDGmI6+65u27GUziQgEixvTlZdTFtsHTgkxIXuk4b0/F4C0RQEaATEpeMfq7XU+JADICZNLyGqlL5EIBIsaE5V0ZukvjQgEixnTlNVWXyIUCRIypyttVF9sGTxAyTXmNKwZviQAyAmSS8lqpS+RCASLGBOW1qRi8JQLICJDJyWtXMXhLBJARIBOTd4C6RC4UIGJMS15WXWwbPHnIlOQdpi6RCwWIGNORdzXUXRoXChAxJiPvYHVrSO+Tx8aZABILMhF5OXVt1x/7avDksXEmTgGIFWQS8g6vGBoI5KUImYK8grqWyz5+7Xvy2DgT1wDECkJfXqeKoYZAXpIQ6vI6Vgx1IpCXJIS4vJeOFUOTCGpeipDQ8gbdAbHdurIMl/0ra3l9JIWIFoR73uP+f0W4bhuMcV56ELrytvv/leGytQSJCwWIGFTlPe7/VwS2DZ4lhKa87RhDxcC2wfOEkJT3ODxWMLBt8FwhBOVlR3a/YtvgGUPIycvflcC2wXOGUJOXvxmMbYNnDaElbwh1iVwoQMSgJO8qjLs0LhQgYhCSV6YuiTYGJBCEjLwSdbFt8MwhRORVVQwk2hiQQBAa8iqLXRJtDEggCAV55RWDHUMXgBCFpC+vdoyBRBsDEgiSurw9w2Mk2hiQQJDE5e2oi22DAWEiaXn71CXSxoAEgiQsb7diwLbBgPCRrrwG6hJpY0ACQVKVt6sutg0GRIg05TWpGPoY5gEIUUiS8hp2u1qGRQBCFJKgvObqEmljQAJBkpNXqBiwbTAgikhNXit1ibQxIIEgackrVRfbBgMij5TktasY5IwhAQhRSELy2qtLpI0BCQRJRl5B3VLePj6JNgYkECQRecWKwSxItDEggSBpyDtQXSJtDEggiJW8D7dZlr24U78mvMFI3sHqEmljQAJBbOR9vMqKWHxQvSa+wUDeoRUDy3ALQIhCLOTd32RP7g4P+Z/38tckb+iX10VdIm0MSCCIhby75cmnQ9m9vpO/JnlDn7w6dU12iSDRxoAEgljIu8me1n+/lL8meYNeXl3FYLa5FIk2BiQQxELedd2hbpmygHtN8gatvNqKAfIC0gcxlzevaCs3d8vWTe412Rs08q4udcWu4Z6UJNoYkEAQc3kfr+pRhLq0FV7rvuFZEar93y75DQCFaOUNsPccYiIRUF7230Y3iopB++8SPS8gvZBR5K2KXf2poeYFpA8yQs3bjDFAXkDcINFHG47DY32nhnFeQPSQ2OO8zPCYj1Mj0caABILEvcPG76HWf8jeINHGgASC2M1tOJPMbThj5zZ038DJ27mhBnkBcYPYzCrbLZlJY+uqRuBe477h8EV0b6hBXkDcIIPn89byGs/nFe8FQ15A3CCRnqSQTcGBvIC4QeLIK52CA3kBcYPEkFcxewzyAuIGCS+vctIu5AXEDRJcXvWkXcgLiBsktLya+eaQFxA3SGh5EYhwEVZerdjxDymPZBJBJpIwyATyphDIRAzIq41kEkEmkkhTXgTCT0BeBNmAvAiyAXkRZAPyIshGYHntV/aNmMjnV1m2qF+rV2nN2of3Y2bCHzxWk4gH2t9kTRSpxGuTQ3Gwp5rs5G0SVt4BK/vGS+RjdWkW5YN39cMg4S+ULBPu4LGaRHKgjrzR2qSI5iEHaXaKNgkq75CVfaMlss0Wbw/Fa+W12XbaLmYm3MFjNYnuQLtlqUmsNimSWWf8wYw0CSrvkOeOYyWSN0j5oH792pp5rj9yJvzBYzWJ5kBN08Rqk7wLuc468hppElTeISs+xEqkXaSqvET7m/D/Tasy4Q8eq0k0B9pUPVy0Nsm7+OziV15eI02CyjtkrZ1oiRx/9rJQ+eSn/F//2YeQWagy4Q4eq0nUB2o6uGhtctieve/WKEaahJR30CpnsRJpouqBm88mobs7aSbswWM1ieZA61qRWG1SBS+vmSYh5R20vmSsRJqo/o/M/+PKPxJ8uw39IV+aCXvwWE2iPlBbWcZqk0N9NFZeM03mLu+2ujR1mScM2ETJhD34+PI22URrkyogr0UiVWyXC+4DbOhKU3vGxcFHl7cZauDTCplJc5S05E2/5t10/0sMrYz2jIuDj17zii0Quk2qSK3mTX20Yf9RKOeCXyjdGZcHH3u0YSMUCWPIO/5oQ9LjvNXilp/aLyMpI8mEP/jY47ztnYl4bVIfJbFx3pTvsBV3JNn/D8rmEQu+GJlwBx/5DltbX0ZskzI68o5/h23Iyr7REtmw3+yW2cV9cZcy9P+Qsky4g8dqEsWBmBohWpuU0ZHXSJOws8oGrOwbK5Fmvl99U32TRZrLJWsS7uCxmkSaCVcjRGuT6sC1vBaaRJzPa7ayb6xEtq27VVoPxeTeiwjFnbRJ2IOPMp9XJm+8NjmI8o4/nxeBCBiQF0E2IC+CbEBeBNmAvAiyAXkRZAPyIsgG5EWQDciLIBuQ12+0N+4Wz9/3v7tY5aN8/FN2D3b/87/0/N7cA/L6Deaus4FeOnk/X6vnlUHeMiCv32Dl7Z/UqJG3mDgEefUBef3Gtl7a69vH7howstBICHn7A/L6jUbe0q/ywYjs5f8si+moh/9lFqUsp0kt3vxfp+d9uF1m2fdv7ou5Ve2EN/XvzT0gr99o5S30y7/K5f0xa74qP8iV3em2mqB6ds3JW79aPbbVyKv5vbkH5PUbkp63iN8xwxAfuJnwjLzHV58e5dX93twD8vqNtua9PfabT+9bl4vlEOtXq0dsWHnzVxfvy3fmiLrm1f7e3APy+g12tKERMv+reJClXvM2KxciOD7r0srbfArb/fjmt/YDm/b35h6Q128w8pZV6iarnqthXi9XFqnc44bK+PGF+jvt7809IK/faGX7vhofkMi7+NBqaiWv+HtzD8jrN7adLRwYeY+v2/W86t+be0Bev6GSt6hS2WVGNDXv4+9f/Bdb86p/b+4Bef2GSt5yDOHuUH/kakYNroTRhrfFCmqFr8xog/L35h6Q12+o5LUb581/pfrmyT3GedUBef2GUt7mdkUlXbMEzA/SO2zlN+vaYs3vzT0gr99Qy1tuuJk9Z+covLhby+Y2VGvU7G/r5WrUvzf3gLwIsgF5EWQD8iLIBuRFkA3IiyAbkBdBNiAvgmxAXgTZgLwIsgF5EWQD8iLIBuRFkI3/B5Mbs2DUxdKXAAAAAElFTkSuQmCC" />

<!-- rnb-plot-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuTkFcbmBgYFxuYGBgIn0= -->

```r
```r
NA

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Calculate hosmer lemeshow statistics

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuXG5obC5vcmcgPC0gcGVyZm9ybWFuY2VfaG9zbWVyKGdsbS5vcmlnLmZpdClcbmhsLm9yZ1xuYGBgXG5gYGAifQ== -->

```r
```r

hl.org <- performance_hosmer(glm.orig.fit)
hl.org

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiIyBIb3NtZXItTGVtZXNob3cgR29vZG5lc3Mtb2YtRml0IFRlc3RcblxuICBDaGktc3F1YXJlZDogNS4zMzNcbiAgICAgICAgICAgZGY6IDggICAgXG4gICAgICBwLXZhbHVlOiAwLjcyMlxuU3VtbWFyeTogbW9kZWwgc2VlbXMgdG8gZml0IHdlbGwuXG4ifQ== -->

Hosmer-Lemeshow Goodness-of-Fit Test

Chi-squared: 5.333 df: 8
p-value: 0.722 Summary: model seems to fit well.




<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Create smoothed cal curve from Van hoorde et al

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucDwtQ3JlYXRlU21vb3RoZWRDYWxDdXJ2ZVBsb3QocHJlZGljdGVkLCBvYnNlcnZlZClcbnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApXG5TYXZlU3RkU3F1YXJlRmlndXJlKHAsICdmaWdzL3Ntb290aGVkLWNhbC1jdXJ2ZS1vcmlnLWRhdGEucG5nJylcbnBcbmBgYFxuYGBgIn0= -->

```r
```r
p<-CreateSmoothedCalCurvePlot(predicted, observed)
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/smoothed-cal-curve-orig-data.png')
p

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArwAAAGwCAMAAAB8TkaXAAAA+VBMVEUAAAAAADoAAE0AAGYAAP8AOjoAOmYAOpAAZpAAZrYyMkkzMzM2Nk06AAA6OgA6Ojo6OmY6OpA6ZmY6ZpA6ZrY6kLY6kNtHR0dNTU1mAABmADpmOgBmOjpmZjpmZmZmZpBmkJBmkLZmkNtmtttmtv+QOgCQOjqQZgCQZjqQZmaQZpCQkGaQkLaQtpCQttuQtv+Q2/+lpfGzs/+2ZgC2Zjq2kDq2kGa2kJC2tpC2tra2ttu227a229u22/+2/9u2///bkDrbkGbbtmbbtpDbtrbb27bb29vb2//b/7bb/9vb///r6+v/tmb/25D/27b/29v//7b//9v///9KLehjAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAgAElEQVR4nO2dbWPcNnaFKVmRa8ezSWuvI61dt2mi1t02qTVyt922rqfJNrv1jmRp/v+PKd8HIAHwXrwRIM/5kGhGM4e44GPoEgRxiwMEZapi7gZAkK0AL5StAC+UrQAvlK0AL5StAC+UrQLB+yUEhVNgeA2/+7MHfx8eMMnS5OLiAvDCJEeTEt2LA+CFSX4mFzW7gBcm+Zk06AJemGRn0qELeGGSmcnFkV3AC5OsTAR0AS9McjKR0OXDe3/5RH7j7m1RFM8/KF6I9lOtslWafQyTICYXA3bZ8N4UMrz3l0Wlk3ejF5K9uVX2SrGPYRLGZIguF96Hm0KG9+GqePzhcFf+99PghWxvbJWDEuxjmAQxGaPLhPfuVTGA93bz6OOhHnK/H7yQ7Q2tclJyfQyTICajjKEWB959Ubz4SYZ3177cFS8HL2R7basclVgfwySIiRpdJrznP5YAS/DetIPsvkoVpBeyvaZVzkqqj2ESxkSDLv+CTYa3zHIbXm83jz9JLwb26la5K6U+hkkQk2906LrCe3/ZzixU6a704tCu5f0zBDnom0q6XwaEV/y3oRJGXphMqRx1v9H/FvDCJFmTOtk1mCDnhUmiJhfPSm23weDFbANMgphst9sG3ZDwYp4XJr5NKmJ7dEPCiztsMPFqsu3RbdkNCG+Z556LaxvOsbYBJrYm204CumHgbZeX3W6EhWTSC8leJcALE0FqdIPCi/W8MPFgspXQFdn1CS9TgBcm09rK7G63gBcmWZhsTeienZ0BXpgkajJE95kE7hnghUmiJnJ+oEQX8MIkSRM9uj25Z8h5YZKeyXjUfTYkFxdsMEnRZDtmdzzmAl6YJGayHUmBrvBbwAuTREzU6D7TkQt4YZKGyRhcYdhVkwt4YTK/iRJcArqAFyazmmjAPWYMenIBL0xmNNGSOxx2dZ8CvDCZyWQSXeOoC3hhMpfJtRndZwR0AS9MYps03GnhldE1sgt4YRLVpONOBy951AW8MIlpInKnhldKdqfQBbwwiWQy4E4Fr5QxTKMLeGESxWTEnQJeLrqAFybhTVTcjeBlZgyAFyYxTJTcDeBtMgbidRrghUkMEy13MrxW6AJemIQ00XMnwsuaHgO8MAlvYubuCK84x8AhF/DCJJDJFHc9vPboAl6YBDAhcHctomvJLuCFiW8TEnfXHbrMqV3AC5OAJjTurrc2dyUAL0zCmZC5u7afYwC8MAlgwuDuqeXUbirwRqw3B0XQNUNPS11fN+hyvifL0BiMvDAhm7BGzHLQfeqU7M4/8hp+B3jzMuGi++zZte30GOCFiT8TLnHOcwyAFyZeTNi8ud2VALww8WZig66vYRfwwsTWxAY2v+gCXphwpX9qnYHuGdlkeHTACxMLdcjYwTt4VoJmMtEUwAsTkxREWcE7vKFGMJluGOCFiUY6qCzgHd8LnjRxDAfwrtNkCis2vHXGMFjGMGXiGg7gXZ0JiUUmvGKySzZxDgfwrsyESCMPXsWoO2XiIxzAux4TDo4ceDXoGk08hAN4V2PCIZcFryrZnTTx0yeAdw0mTHI58B7RZZh46hPAu3gTPrl0eEtyT/V3gjUm3voE8C7axApcA3djdA3sakz89QngXa6JNbla7mRNoKsx8dgngHepJtbc6rmTVJJ7OrF0TGXis08A7zJNrKk1cCfq9LRll2nitU8A7wJN7Jk1cSfojIKuysRvnwDe5Zm4QKvnrtdZhy7fxHOfAN6lmVivIzdz16lB95T0kMTQxHefAN5lmZi440hjctYOu1Ym3vsE8C7HxMwdTwqTemqhQtfWxHufAN7FmBi442tk0qP7zNbEf58A3qWY6LmzkWTSLrs546E7MAnQJ4B3GSY6ZGzVm5z1Om22ErEx2dqzC3gXbqJHxlqNyZHcs3YrEQuTWkH6BPDmb6JHxl7XErn9Dk5sk1Zh+gTw5m5iQMZBIriW6AotCdQngDdrEyMy1pLJ3dpkDHJLAvUJC967t0VRPP9wfOPhquj06OPhcH95/Fm2Z7aKrIy582AygYydBuRaD7vHlgTrEw68LZsn7/p3BvDebgBvRBMzMjY6G5HrgG5S8JakPv5wuCv/+2n4q9tNTfS+eKKxZ7aKrFy582EygQxbErlC/T9rdluTcH3CgPd2U4+o5fj7/eA3JdYvq//fNP9T2DNbRVam3HkwmUKGqeGY25g4oduYhOwTBry7dlzdjRDdNYPxw5WQUcj2zFaRlSV3HkymkOFpnCwI9f8clBC8N+2Iux/mDd1YfH/56IdXRXEuEAx4g5hMIsOQgtxtV//Pjd2qJUH7hA5vmRs08N5uBvDetDR312vCyAx4A5hMI0OVGtxKT93RrVoStk/o8N5ftklBm/sKv+iH5OqK7vPbdj7iy0qxitutSPb1+GQJae7od039P2cF7gkP8O66NKL74eY46YCR17cJZbyjSDvmbrv6f+4K3Sfu8HZTDUcJSTHg9WpCQoYAr4ncrv6fu4L3iXvOOxyIpXcAr0cTIjOT3JnI7eYYPMAbvk/cZxt2ozsTgDeECRmaiafWKeh6gDdGn7jP8/Z3JvqhGWmDfxMGNibuzOQK02NLg1d9h61PhfsLNTEJBrx+TDjY6LmbQFe8K+EKb4Q+4a5tOB+vbRByhNtN8eLT4e4VFub4NmFxo39q3UiufEPNEd4IfXLgrSprb0I0I203HyZmwLtiuOwM8LqbcMFRcTeR6W5HS3Cc4A3fJ4048ErreVXwHu5el+i+EAZmwOtswkZnzN00uqN1DC7whu+TVix4+QK8riZ8dgbcEchVLMFxgDd8n3QCvEmb2MCj3HLB9AXVEhxreMP3yVGAN2ETO3wE7ijoqlc+WsIbvk9EAd5UTezoEbgjkatbtGsHb+g+GQjwJmpiBY/AHR1d5cpHG3iD98lQgDdJEwt0RO4oF2mV9It22fCG75OxAG+CJlxwBtwRyTU+5sOFN3ifqAR4kzNhYjMSA139sxIceCP0iVqANzETBjRK0cidfLqSDG+MPtEJ8KZlQmVGLeKgS3gwmAJvrD7RCvCmZOK0ooB6lUbaSsTckph9AnizMHFbUdCjO2lCeTBYbRK/TwBvHiZaZEg6DroTJrRn2kMXr/RiAngTMdEhQ5KULhhNqFuJhC7k48UE8CZhQuFOq0GqazIhbyUSekd+LyaANwETFTJ0DS/T9CaMXXA6k9n6hGICeOc2USFDl2KCQWfC2nysNpmtT4gmgHdeEwUyDCnnxtQmzH3zrmfsE7LJGN7Pv8j6o8uRAa9RI2QY4GoXjilNmPvm5dGxI3i7uhLSfv3WArwmjZhhwau7I6Ew4W75mEnHAt65TFTQsHcnJZlwd9qdr0+YJpPwfvE14A1iQuJOK8N94KGJBbqZdKwm5/2H4uS7Px4eftooqqdwBHg1onGnk3ENg2xihW4mHauebdgXx/2ihwUoWAK8apG408q8/EY0sckY5uoTvokS3oerLtMtkwinoRfwqkTiTquppWOCiS26mXSsEt4S2SO8uGDzbULiTqfpVY+9CX+OYb4+sTLRwdtmCzvMNvg2IXGnFWHFrmX9vzn7xNJEnfPeVHuS/f73//66/P9w62iWAO9AZoAm4KUtNreq/zdnn1ibqOHti1JJWz5aCPBKmkKIsqn5JIcW9f/m7BMHE83ahp82XtgFvJJI3OlEfj6NX/9vzj5xMdEtzPn8u19Vdyi+dZrlBbyiCBQZ4CWSW5m4XKhF7hM3k9CrygKXkctHbtX4mlGX9FFu/b+5O8ZBeng//+E/HIfdA0beXrQxUDfy0ofdqv6fY8YQrU+cTXTwVknvo4/3v/k7tyMD3kZElNTwUpPdbXOhlsVmN15MNPC+b9aT3V+6zZQB3kZUlvTlJEhf59b/m7VPPJho1zZ08Ep1q9gCvJXINCm4o7PbzTFksVOTFxPtTYonf7p89LG6w4a1DY4mVJaU3HGHXaWJSvP2iR8T3e3hk3f3FbxY2+BsQiZXwZ0FuplsM+bFRL8wB/B6MWGgO+SOnzEQ4Z27T3yZTIy8e6QNbiYsdlWFfKjsakzs0M2gYytpc97H/3n56L+rCTMszLE34aGrKOTDR3cS3rn7xKMJFuYENOGyK3BnlzEQ4J27T3ya6Bbm9PA6PQW0anjZ6ArcWQ+7WyO8s/eJXxPdHbaqinCp8x/djgx4reB1QNcAbwJ94tfEsDDnl19+cT7yiuG1YLfjjsiuZuWjBt4U+sSziXq24bcffBz2sGZ4bdgV6v9Zo6uBN4k+8W2imypzXskr26u0aHit2C25Iy85Z9X/S6NPvJvo4PWQ74r2zFaRlWof27Hb1/9zQFcBbyJ94t9ECe/D37dTDc9d04eVwmvJLnVyl1X/z0M4tprpgu3h51ftNK9b+rBOeF3YJXyMUf/PSzjJmuhnGx7qp9iwbwPfJCy79Pp/nsJJ18QwVXY4fH4LePkmQdml1//zFU7CJoaRt80cAC/LxBbdml1/9f98hZO0iQ7en18j57UxcUF3snglsf6fx3DSNlHPNrzHbIOdiRu7Xur/Jdcn4Uwwz+vTxJFdD/X/vIaTuokO3pNvnaoADe1VWh68LuzWP7jV//MdTvImGnh/7eOwB8DLY9ep/p//cJI3Uee8V46lKEb2Ki0OXmd2dfBy0E2sT4Ka6NIGpwmysb1KS4PXHl1zzetpdIOEk4GJbuR13Np0aK/SwuD1wK6u/p+R3UDh5GCinufdb4pzL0t61wOvVcH24V21sckEusHCycJEPfL+Y3Nz7eSrUigiSJFFwXbFk2q8+n8Bw8nDxDDPi/KtdPELtqvYVdT/4w67fsLJxATw+jAZcWfHrmxih24yfRLBRJ02/OH3RwlbTN+9Hd0y7jhvEB99YB3wjrijsjt8j1z/L2g42ZjoFuao1IIqTUR025PU8I4/sAp4h9xZsyvX/9OzGzacfEz08I629a9uXXw43Mk3MPbCdlCKD6wB3gF3ZHQVi3c7Ewd00+iTOCY6eBXb+t9uuuFV2G/6RthSR/EBwMthl1L/L3g4GZlo4FVt679rf94JvD5cCTmC4gMrgHcwaLqw29f/07EbIZycTDQ3KVTb+t+0P+6FtOD+8tEPr4ri/J3mA8uHVx40yeyqf3NtHnYjhJOViRpe1bb+ZUbbsHm7OcLbbyf5Uv2BxcMrDZqu7Jrr/8UIJy8T3TzveFv/+8s2Q2hT21rV5tMfqgc1y98NP/BlpbnrzAUWp1pfK0M1wKemAoBzh5qg9KvKSPDu2oH5pkx3VR9Y+Mgr/8V3HHeN9f+ihJObycTIK2zrr4S3U5Xmrg9e6S8+GV09u/r6f1HCyc5Em/OOtvVX5rydKl5Xl/NKeJHg1bNrrP8XJ5z8TNTwKrf1V802HL9QDrYrm22QAaPAa0RXW/8vUjgZmmjmeVXb+iumcfvBtuZ1XfO8A8QI8E4MuxqTSOHkaKK7w6bY1l91h+2m4bWE+OXa7rCx4dWxa67/FyucHE0MC3OG2/qXhJ4P1zaU+cWLT4e7V/WchOIDy4V3CNkkvBp2zfX/ooWTpYl5VdlnKbltM+GT9nZaPebuxMxY+oBkz2wVWbP18QjCCXh1l2rm+n/RwsnTRJs2/PZj9SSb/CibuFy3hbfOL05efBp/QLJntoqsufp4TKEZXg27E/X/ooWTqYkG3l2ZBzTjqNsz8AuFV4HnNLxKdE31/6KFk6uJdmHOo483HqoILhNeFZ5GeEnD7sAkXjjZmmhvUvzlny6L4vEv/4rC2SMp+TTBqxp3p+r/xQsnXxPd7eHv66uv75X3ghlaIrxqQA3wKtidrP8XL5yMTTTwNqsh66VigFeWhlAtvKprtcn6f/HCydlEC+/DVb0mByPvQDpGdfAq2J2u/xcvnKxNlPCW4H5XZQ0vDw83yHkl6djVwTtm1/SYz3pqoXgx0V6w1TcbqseAMNsgigkvZ9htTaKGk7eJGt5mA4Yn7nudLg1eLXZqeJXDrt6jMokaTuYmmpsUd6+KatHC/eW52z6964Z3yO7kdqXX27jhZG5iXtvw8IvybbqWBa+ZOxK7Jos1VaHyYmKG11lrhnfALqUkShbIpGOihffun968eeNcEmhR8JoHTTO7pCJqeSCTjokh523KCLoVVlkSvGbwhvCO2Z1Cd5sJMumYmGYbmmcpMNvQaIK8AbwSu6j/F8bEMM/79ZvXWFXWawo9GV6RXWrZ1ajhLMJEd4et+IsqX3i4KXCHrdYkexK8Q3an0d3GDWcZJpq1DUW/fwi29a80DZ8Ir8AuEd1t3HAWYqIbeVtkS4wx8h6Y8B7ZJWYMKywb7MVEnfPuupF3Xwgb9FpoIfAS8DvC27NLRncbOZylmKjhbR5iPxzuLh1LYS4DXgp/PbwSuyR0t5HDWYzJCN7711XlwE1dQrD679d/s/q0gQRgB2/HLh3dbeRwlmMyhleswYY6bAdqMfdriV16xrDSgu1eTADvlAkRweshu8TvbSOHsySTEbxSAcFBEUG+8oeXyuD1kV0OutvI4SzKRLO2wZeyh5cM4XX3yA8nYxiuPc8CmXRMAK/RhAxhBy9v2N1GDmdhJhp47/7tzZs3//xh9HGu1gRvxS4P3W3scBZmooT351+1F2tffOd45MzhZXBYVfk5dWQ3D2TSMVHB+16Ya3jhduT1wNuyy/jGNno4SzNRwLuTZsqcVkRmDi+L3dNTV3QzQSYdkzG81W4jJ99W82Off7cpXG8Pz11nzkWcyoCnp6f6+n8oChhGY3h3ArDNtjku8Bp+l/rIyxl3S3ZZw+42fjgLNBnfpLiSKwCtdkkkg8QS3VNO1XfA68dEdXtYyBTWuxidzuGzil1O1Xctu3kgk46JCl4BV9f9nrKFl8fuGafqu57dPJBJxwTwKk3o6Jbs1mtxOPBGD2ehJsqc91gkcL/OnJeObscuB97o4SzVZDzbID4wXK2PXOFjQBx2u+cm6PBGD2exJmN4y8G2eNysaqj3zfle8S2ylgxvdS/4+KwlGd744SzWZAxvlTcUxclfvfnrTbNJr4uyhJeMrsAuGd744SzXRHF7+P6VcHd4hfvzktmV9hYhwhs/nAWbqBbmPLzt2XXcZy9HeMnoyvuR0eCNH86STVTwVvubfrUpTr76W+cFvfnBS0R3yC4N3vjhLNpEDa83ZQcvkd3taO9oCrzxw1m2CeCVTMjojvbsJ8AbP5yFmwBe0YSErpJdArzxw1m6CeAVTMjoKuqrTcE7QziLNwG8RxMKu1sNu4B3BhPA25uQ0VWxOwVv/HBWYAJ4O5NpdE3sTsAbP5w1mADe1mSa3e5nJbtmeOOHswoTwNuYkNHVsGuEN3446zABvLXJFLqT7BrgnSGclZgA3spkil0JXSW7enhnCGctJoD3cDD9xR9s4KRlVwtv/HDWYwJ4TX/xh9uV6tnVmUQPZ00mgNcArwJdHbsak/jhrMlk9fDqB83hlo9GdtUm8cNZlcna4dVyN9rg3Myu0iR+OOsyWTm8Wu5GO+2a0VWZzBDOykzWDa+Ou/Em0VPsAt4ZTFYNr4Y7RUmUSXbH8MYPZ3UmK4ZXx51ib/5pdkfwRg9nhSbrhVfDnaqsBIHdIbzRw1mjCQveu+qZ+OfyI8U/vy6Kk/a9rnrmcW++ZOHVcKcsokZhdwBv7HDWacKBt2VT2ui/Lb5yUm8KdbvJBF4dd8pqPiR2ZXgjh7NWEwa8D1fVHmZ3V+K+kfvi5LtD9V7N6360OVSa8Gq4UxeiorErwhs7nNWaMOC93dSEluNvv/deyXNdAqB972ZUwCJJeNXcacquEtkV4I0dznpNGPDu2nF1d0T0/rLNIWpsH65GpYMShFfDnab+H5XdHt7Y4azZhAHvTTvi7hX7Tdfw3l8++uFVUZwLBKcHr5o7A7o0djt4Y4ezahM6vGWK0MB7uxnB24zA3fWakDwkB6+SO22hdjq7Lbyxw1m3CR3ePkVoc19Ru3ow3te7Un9+285HfFlp7jpzstTl/J4+1dT/q9hFXcB05QXefYPrrs0nbo6TDmmNvMoh89mzp5r6f4xxtx55Y4ezehMf8O43J9Le/0JSnBK8GnRLeD2wW8IbORyY+Mh5d8P6xALdCcGrQ/eZej0v41qtgddDNHkgk46J82zDw/tRbe0U4VUC116nqeDlsrvN4mwvzMRtnrcej7uyFf3QnF7aYEJXCS+X3UzO9sJM3O6wHR5uhGG4vVDrbruJ9sxWkUXz0KHbXqiN4GUPu56igQnPhLe24Xy4tmEnvrjdFC8+VcXb0lqYo8RNnNkdwmvDbh5ne2EmnFVl7U2IJsWth9luDWRbr21XDJedzQ/vJLojePkpg69oYMIz4cArreet4d0XEryHu2px7wvhcm5ueLXoilO7Mrw2w66naGDCM2HBy9e88KppG90LluC1ZDePs70wk+XCq4FNsYxBhNcqZfAVDUx4JkuF14Du6GawAK81u3mc7YWZLBNeHWzq1WNHeHnseo8GJjyTBcKrhU238LGH14HdPM72wkyWB68JXfXysQ5eFrshooEJz2Rh8Oph06Lbwcub3g0SDUx4JouC1wrdFl4Wu4GigQnPZDnwGljTZwwdvI7DrqdoYMIzWQq8JtbM6Fbwug67nqKBCc9kCfCaWZtCt4TXA7t5nO2FmeQO7xRqExlDJR/oZnK2F2aSObxTtdan0fXEbh5ne2EmOcO7nSpXTUCXxa6pMVmc7YWZZAtvg9NU/T8Su1PDN4XdPM72wkyyhPfIk4E7+rBLg3eiTVmc7YWZZAevTBS1/p+JXRq8U+3K4mwvzCQreMdIGer/kdmlwDvdtizO9sJMsoFXDZWaOxa6FHgJ7cvibC/MJAd4DVSR6v9NsDsJL6mRWZzthZmkD6+RK3X9P+Kw282QTcBLa2YWZ3thJknDO8kgpf7fFLoT8FIjzeJsL8wkUXinsVVyx8gYhBsTBnjpkWZxthdmkiK8VHKH3Nmha4KXEWkWZ3thJqnBywB3wJ1VxmCGlxNpFmd7YSapwMuEdsyd7bCrh5cXaRZne2EmScBrSa7AnQO6OniZkWZxthdmMju89uAeuaNlDNoFZCp42ZFmcbYXZjInvE7YCtw5DbtbBbw2kWZxthdmkj+8ruiO4bWKNIuzvTCT3OHlZAy6NecyvJaRZnG2F2aSObxP3dGV4bWONIuzvTCT0PCaChhyiksqpS1dKalBl+YYqXQj5EUZj7yG+n+MYVcceV0izWKoWphJvvDq6//x0EXB9mxNcoXXUP+PiS4Ktmdrkie8/RyDCd6zMxq7tYlzpFmc7YWZZAmvof4fH13UvM7WJEN4TfX/LNBF2eBsTbKDV74roYaXQ+42kxMFk7Fyg9dY/09il2h4yOREwWSsvOA11/8TyGXtPpbFiYLJWDnBq1jHoC7Yztw4L4sTBZOxMoJ3ov7fEV1ysts1xEekMIlvkg28U/X/uAnDVmiIj0hhEt8kE3h1Kx/l4pVkdgcN8REpTOKb5AEvqf6fJbqZnCiYjJUDvOb6f2dnLHK3iob4iBQm8U3Sh9f4rMSZO7qZnCiYjJU6vOHRzeREwWSsxOE1oMvkVotuJicKJmMlDa8eXT65W0NDfEQKk/gmCcOryxjEZIFayGdramQWJwomY6ULrxLdM5HcLbEWytRS8yxOFEzGShXeCXS7d1DIZ80macKrzBhUaS5qoazZJEl4x+jqLtBQC2XNJgnCK6F7Jmn0WezIv2aT5OA9ZgxnZxPkmuDlNDKLEwWTsVKDt0OXAK4BXl4jszhRMBkrLXgrdGnU6uHlNzKLEwWTsVKC97QSb70CttZds0kS8DasyugSyJXhtW9kFicKJmPNDm8Ha48ujVoZXrdGZnGiYDLWnPAKGUKdMbCwbeH10cgsThRMxkoDXtre/EMdMuljmAQyYcF797Yoiucf9O+NPjAN75ZaEmVMrjkyumCSqQkH3vvLotLJO9174w8QLtiIJVFkaAmR0QWTTE0Y8D5cFY8/HO7K/35Sv6f4wDS8VHTZkdEFk0xNGPDebh59PNTD6/fq9xQfmILXjO50m7LoY5gEMmHAuyuetP9/qX5P8QEzvIqMYboZsgfz8zBZkgkD3pt2QN0LaYH0nuIDRngvSk0f1qgs+hgmgUzo8JYZbcPm7aZnU3pP9QEDvBffOKObSR/DJJAJHd77y3YWoU1tR+8NP/BlJV39t28qRaw3By1QAeEV/20MVWUMPv5dZjFAwCSQySzwNsku4IWJm8kMOe9Fe6EGeGHiZhJ9tuGin2MAvDBxM4k9zytMjwFemLiZxL3DJs3sAl6YuJnw1jacK9Y2nItrG4YfkOC9kO9KAF6YuJkw4C1HVmHR2E2TI0jvSS8k+0rDG2qAFyZuJhx4peW6Lbzk9bzje8GAFyZuJix4+ersLxTrGAAvTNxM4sCrXIIDeGHiZhIDXs3qMcALEzeT8PCqMoZagBcmbibB4dUv2gW8MHEzCQ2vYb054IWJm0loeCEonMLCawQ7/iHVSqYhaIlChJYA3hSElowFeI1KpiFoiUJpwgtBfgR4oWwFeKFsBXihbAV4oWwVGF7+zr4RG/Lz66I4ad9rd2kt+of3Y7ZEPnisLhkf6OGq6FQ1JV6fHKqDPTG0Tt0nYeG12Nk3XkPeN6fmpH7wrn0YJPyJUrVEOnisLlEcaABvtD6p1D3koGydpk+Cwmuzs2+0huyLk+8O1Xv1udkP+i5mS6SDx+oS04FuNzUmsfqkasxNIR+MhElQeG2eO47VkLJD6gf12/duhOf6I7dEPnisLjEcqOuaWH1SDiGvigG8JEyCwmuz40OshvSbVNWn6OEq/J9pXUvkg8fqEsOBds0IF61PyiG+ePGTDC8Jk6Dw2uy1E60hx9+9rFB+9EP5r/88+NlStUQ6eKwu0R+oG+Ci9clhf/7jMEchYRISXqtdzmI1pFMzAnfXJqGHO2VLxIPH6hLDgW5aRGL1SSMZXhomIeG12iVd53IAAAOESURBVF8yVkM6NX8jyz9c5SXB57ehL/KVLREPHqtL9AfqM8tYfXJojybCS8Nk7fDum1PTpnmjCZsoLREPPj+8XWui9UkjwMtoSKP95kS6gA2daRojrg4+O7zdVIPcrJAt6Y6SFrzp57y74Z/E0MgYI64OPnvOO+6B0H3SKLWcN/XZhof3o3Qu+IkyRVwffO7Zht0oSZgD3vlnG5Ke5202t/zY/xgJGUVL5IPPPc/b35mI1yftURKb5035Dlt1R1L8e1B3zzjhi9ES6eAz32Hr88uIfVJrAO/8d9hsdvaN1pCd+OJ2U7z4VN2lDP0XUtUS6eCxukRzICFHiNYntQbwkjAJu6rMYmffWA3p1vu1N9V3RaS1XKoukQ4eq0uULZFyhGh90hy4hZeBScT1vLSdfWM1ZN+z2zTrrlrc+yJCcqfsEvHgs6znVcEbr08OY3jnX88LQQEFeKFsBXihbAV4oWwFeKFsBXihbAV4oWwFeKFsBXihbAV4/aq/cXfy9Y/Tn652+agf/1Tdg3343b9MfG/tArx+Jdx1JuBlgvfnV/p1ZYC3FuD1KxHe6UWNBnirhUOA1yzA61f7dmuvz++He8CoZIAQ8E4L8PpVB2/NV/1gRPHyfzbVctTD/wqbUtbLpE6+/b/ByHv3dlMUX3z7qVpb1S94039v7QK8ftXDW+FX/lTC+1XR/VRfyNXD6b5ZoHr+SoK3fbd5bKuD1/C9tQvw+pVi5K30a2Ea4p20El6A9/jukyO8pu+tXYDXr/qc9+1x3HzyqWe52g6xfbd5xEaEt3z35Mf6k6VFm/Mav7d2AV6/EmcbOiCrjOB202YGFYsVkf2zLj283VXY7Vff/rG/YDN+b+0CvH4lwFtnqbuiea5GeL/eWaRhT5oqk+cX2lfG761dgNeveti+aOYHFPCevOsxZcE7/t7aBXj9aj8o4SDAe3yfN/Lqv7d2AV6/0sFbZaniNiOGnPf+N8//S8x59d9buwCvX+ngrecQPhzaS65u1uByNNvwXbWDWsWrMNug/d7aBXj9Sgcvb563/Erz4vEnzPPqBXj9Sgtvd7uiga7bAuZXyjts9YublmLD99YuwOtXenjrgpvF1+IahecfblRrG5o9ah7ettvV6L+3dgFeKFsBXihbAV4oWwFeKFsBXihbAV4oWwFeKFsBXihbAV4oWwFeKFsBXihbAV4oW/0/nt0TFl7BSv0AAAAASUVORK5CYII=" />

<!-- rnb-plot-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Boot strap performance in the original dataset

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArwAAAGwCAMAAAB8TkaXAAAA/FBMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAZpAAZrYAv8Q6AAA6ADo6AGY6OgA6Ojo6OmY6OpA6ZmY6ZpA6ZrY6kJA6kLY6kNtmAABmADpmAGZmOgBmOjpmOmZmOpBmZjpmZmZmkJBmkLZmkNtmtrZmtttmtv98rgCQOgCQOjqQOmaQZgCQZjqQZmaQkDqQkGaQkLaQttuQ27aQ2/+2ZgC2Zjq2kDq2kGa2tra2ttu225C227a229u22/+2/9u2///HfP/bkDrbkGbbtmbbtpDbtrbb27bb29vb2//b/7bb/9vb///r6+v4dm3/tmb/25D/27b/29v//7b//9v///+8f07FAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAW80lEQVR4nO2dC3vbxplGIcnSSmxcp/aGTJxNs27Mqq03Wa8Ze7uNV2btxGmUkhKJ//9firngRoAkMAA4+IhznsemSAHzQuDhcDCYAYIQQCiB7w0AcAV5QSzIC2JBXhAL8oJYkBfE0kzeFtX/tb2iCBcS3hTkJVwsyEu4WJCXcLEgL+FiQV7CxYK8hIsFeQkXC/ISLhbkJVwsyEu4WJCXcLEgL+FiQV7CxYK8gw5/9+6dx/SmIO+Qw9+9E20v8g45HHnbYaj+IK87yDvocNHuIi/hckFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvL7Dr66uPKYjbxsMVd6rK6/2Im8bIK8XkLcNkNcLyNsGQ5WXNq87yEu4WJCXcLEgL+FiQV7CxYK8gw7nnhStMFh/fIZzN6B2GKo/yOsO8g45HHnbYaj+0OZ1B3kJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5B12uOSzw8g78HDkbYPh+oO8riDvsMORtw2G6w/yuoK8ww5H3jYYrj/I6wryDjscedtguP4gryvIO+xw5G2D4fqDvK7U1W81uWyw9g6G6w/yulJXv1mAvMcUPiB517MAeY8qfDjy3j0NupL37OystbLqM1h5/e72ptTRbxEET953Iu+ZoaXS6jNQeX3v9qbUkvf8ZSQw8h5NuO/d3pS6+nUi79mZ5904THm97/amuMsbKH5thWQvtlMcVKPl3d62mvuh5lVQ84qkF/J6b3wNU17vu70pyKtAXpH0Q17fHY4Dldf3bm9KX+Qdrj+hz7u3DugMW4i8HYC8riCv93DkdQV5vYcjryvI6z0ceV1hJoX3cOR1BXm9hyOvK72R95/tFVUf5BUJ8iqQVyTIq0BekSCvAnlFgrwK5BUJ8iqQVyTIq0BekSCvAnlF0p28/+yWRtu9CfKKpEN5G5W8D+RtB+QthZq3IsjrCvIqkFckNBsUyCsS5FUgr0iQV4G8IqHNq0BekXCSQoG8IkFeBfKKBHkVyCsS5FUgr0iQV4G8IkFeBfKKBHkVyCsS5FUgr0iQV4G8IumLvG2fM6sH8oqkJ/K2f8a3FsgrEuRVIK9IkFeBvJ2xHJ3eRA+rycVt62X3RF7avJ44gLwnr/TDEcvLJU79cAB5vx5HD/OvkbcjkLczlqNvvroN13/6LpJ3fR0EqhGxHAVBMLaP03D5m6hujv5bfvJt9Ot4oQogr/fwduV91y21t2c5ev7nm3D56fuL2/X1ZVQFX9yuJtPo8fTGSDuaJvKOot/HC1UpG3m9h7csb6ulNS89knM+DRfjxcXtQlWokbn3ykxd05oKNpU3kjpeqErZyOs9/OjlXVyGs2kk7zzQjNVteYIgOo6bBfr2PIm86jFdaD/I6z386OVdffWPL26UvLYxsJpE4mpTox+DuPkQy1vjwA55vYcfvbzh999dhqrZoDvNQvVj9J99ErUQcvLGr1cBeb2HH7+886gVsFAHbJG1kZzKT9X9q5u3kbGryTg6Tjsx8sYLVSkbeb2HH7+8ysqF7SpTWkZt3ZP/iY7JbNNXd5l987mRN1moAsjrPfy45e0S5PUejryuIK/3cOR1BXm9hyOvK8jrPRx5XUFe7+HI6wryeg9HXleQ13s48rqCvN7DkdcV5PUejryuFPW7//Gvlcf1IG8LDEbe1mcpbur3fqQmYaw+/4PT2g04InmvumV3+IDl/V4NBI7knehRwnXXbsIxydtqaTVLH668iyCRN6gyEQN5y0DeHN3dDSevn5qX8Y/J6Y2ajFFlRDvyloG8OYqqdiJvVOGevFopefX0jJprNwN52yl9yPKqJgPyNgN5c+yWV7VU4y/5+VhfycGoNw9s09Uusfy00BTYUvMuaDa4g7w5dso7V6LOjGuRnmY2vJqEqf6p+UHpEvPCjOJCm/fib5PT/1cdZlW6G5C3DOTNsUve1edm5o+uYmfT0Mgb/b/6Iv4hWcK8lCWvn6m1NZWmESFvGcibY1dvwyLz/a7kTGpe84vI2cwSs80OsA393ifyHuCqDzmQt53S+3a5p9017yLz/a4atab2VPOHzS9m08wShXbDpn53X2p1z19W2i7kLQN5c+ySN3sUpupYVfPqee/G2ajmzSyxV96Ijx8/Vt0u5C2DZkOOqm1eXfOqZoM6RCtp81aStwbIWwby5qja25C0ec2FInWHbba3YW+btybIWwby5qjcz5v2Nozsxfiy/bz7ehvUkIYETlK4grw5Ko9tKDkNkWVPPy/ytgLylrP3rHBRzwwVzrAhb3OQt5xuh0Suf/xfw4sgeFJ77WYgbzulD1belDnjed1B3gOxRb+oAcHAHFeQ90Bsl5c2ryvIeyDK9VvPOGBzB3kPxPbeBpoNriBvOX9saytitsvLAZsryFvOweQ9qXThBuQtA3nL6VbepJ/3rxXHlSFvGcib449FWtoaBua0H84Vc3IUVUXeNkHeiiBvKcckb6ul1Sxdlrxm2o/pG4jvo12dRL9cTwMDc5qAvDl2yvvJjb1xtgvIG+pDijaLYw5bjr3yrv/0avno6en/fRJf7EY9q1ILI689HG6xPK/y1uQAp8x29TaYmvfTW32T1+jn2TicX+pbvlYh1e/+Y56fa63dGOT1Uvoh5N3xSjzVXYlrZ1uuvrip2vrlgA15O2Z/s+F6Gss70Ze7aS7vfaO1a4O8Xkrvg7zhLJHXzLF0lffHN18bKjWZkbcM5M2xV96osWDlVW3excWtm7zr6yEesB1Vb0NNfMtr+3ljefV1Sh1r3vkgexvaHjKCvDkONLZBVbwnD4PgwW+D4LxWR1tzkNdL6YccXd75kMip+W89kz2et+TT3ia7w5G3nK7l1deFUpc3XY5Ez6RA3oock7zq6pLzQF9RXXSbF3kr0q9JaTXZbPOev1QXNnu+fi37gA15K3I88qreBnMHwYFNwOSATSSF27dGFe5My1vluv7IWwbyHohN/X766iZcvwgGdq0y5D0Ere+XLRcd+VihzbB1bSeQ10vpxyHv6j9eVhO2dO3mMJPCS+lHIm90mPao2j2AStZuDvJ6Kf0A8taeQ1qZzZkUj39wWrs5yOul9EPIW+EVN4rTgE4ev62/dnOOSV6fU99rchzyhuHdd/HNW0+eVZkDFCJvK3RbUe/Bs7zp4Mf9wyCLS+T1S/198KzK4RvytgDy5n/axj55I+7fjIY3nhd5O2SPvMtP/xIEUz3nPTPx3YxKNy+svvg2elJyvfMy/Yy/JfLeqbMX2SYx8rbA0cu7o9mu5B0lU3/Sie/qdq3Ri+aF1eTidnFaMqe4TL/1m4el8tpjupP0+ibI2wJHL++OV5S8ZgZQfuK7EdW+sJpM9YVJmjQb1tfBxdvw7jozZAd5WwB5Y3nTie/mjoH2BX3vbDPDOM+2AzY1NnKD5Uj7rGdalK7dCJ/+nJ2dDTS8X/KmE99tzWteUPLuq3nvXmTMLelrmAeX9jEZcHYU8p4ZBhge9kvezMR31eaNHs0Lq8ll+YT4kmuVlZobqulBpsZdpO0G5BUdHvZI3vV13MVgruUQ9zao/z//vXpQS+TL2ZR3e/9u1OS111EdHZW8Z2ceBfIarpAwtkG3ecvIybvzzISZ4BYmbV/t+q/ySfwZXLjiXXtF7ZXQ9ci0grxf7jmntilvbu3mUPN64ZC7/UCD0cs4VnmH3eYdiLzH2uZFXrHU0O9YexsG3c87FHmPtZ/Xd7jXM2xDkfd4z7Ahr0xq6Keup8PYhvYZjLyt9ynX0c9cCZhRZS2DvK7U0o/xvF2AvK5wNyDv4ccubzu3IiwDeb2HH728FV5xA3m9hyOvK8jrPXzQ8qrxYKc3le//kwd5vYcPWV41PS2cV7/zWh7k9R4+ZHmzNxHMznMPw5KZ7gWQ13v40cu7o7dhfW3OeKWTgOw8dzsNfjfI6z386OXd+coisFN/1GzLqA6289ztrPfdJSOv9/Bhyxux/M0rO5fNWhvOpnbW++6Skdd7+JDlXeiWgb4mQ1zzmnnue2tdBfJ6Dx+yvLq3wR6w5ea522e7S0Ze7+FDllf389rp7ra3wcxzt7Ped4O83sOPXt5aYxu2ThUuAXm9hx+7vAmVzgojr6hw5HUFeb2HI68ryOs9fDDytg7yeg9HXleQ13s48rqCvN7DkdcV5PUejryuIK/3cOR1BXm9hyOvK8jrPRx5XUFe7+HI6wryeg9HXleQ13s48rqCvN7DkdcV5PUejryuIK/3cOR1BXm9hyOvK8jrPRx5XUFe7+HI6wryeg9HXleQ13s48rqCvN7DkdcV5PUejryuIK/3cOR1BXm9hyOvK8jrPRx5XUFe7+HI6wryeg9HXleQ13s48rqCvN7DkdcV5PUejryuIK/3cOR1BXm9hyOvK8jrPRx5XUFe7+HI6wryeg9HXleQ13s48rqCvN7DkdcV5PUejryuIK/3cOR1BXm9hyOvK8hLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFewsWCvISLBXkJFwvyEi4W5CVcLMhLuFiQl3CxIC/hYkFe3+FXVz5vPoy8bTBUea+uvNqLvG2AvF5A3jZAXi8gbxsMVV7avO4gL+FiQV7CxYK8hIsFeQkXS139VpPLBmvvYLhv4XDDm1JXv1mAvIT3hHr6rWcB8hLeF2rpd/c0QF7Ce0Md/RZB8OQ98hLeF2rJe/4yEhh5Ce8JdfVDXsJ7g7u8geJXAEvbau6HmpdwsSAv4WLZr99Ctw+CafwMeQnvCchLuFhoNhAuFuQlXCzIS7hYkJdwsSCv73AmYDrDTArP4Ux9dwd5PYcjrzvI6zkced1BXt/htHmdQV7CxYK8hIsFeQkXC/ISLhbkJVwsyEu4WJCXcLEgL+FiQV7CxYK8hIsFeQkXC/ISLhbkJVwsDeUFiGnJyDr6HT6yHK8bQrhIerPxw30LhxvelN5s/HDfwuGGN0X0xsOwQV4QC/KCWJAXxIK8B2N924cijomO5f3pxSgITh6/3Xx9dnoTLkdx9/aDJ8mbcv/moXrh2c9byljEFwp2YjWJE5/lNJjnLmLVmOUo+vM2+fC0yZY3KWL94Uv1J9s9uJqcvGq6IT2hU3n1TQc14/wvtIKpvEFg3+z198krsc+bZcxKvKhMIm8QXGTtPYC86+tGH7sGRaySPfhEP0XeKkSqnL+MHLl7sWHvaqLUSd/iDyMjT/TuBI9/CFVlG9tVKGM5auDZamIj48Ru6JW80Vrnap/evzZ7EHmrEO21uHqbB7l3c67fhMxbbH+cBcl+XRq7SsqYNdj3ibxR3d+kBt9Dr+RNN8b8zchbhYwf0W7PVL2m4s2+xWaHLkeZ3brQIpeU0aTqTeU1Py2C8WIUnL+yzYY71bp+9NbEm984xZi/bDm6+OVNVOD5y1B99BTTbSEfoq/2B8/N2rkFTBG6uZoWUYtF8vE3ezmWV32XlcXkt6DXdCjvLCPsfb6FaS1MxDSO5lueevWSMtbX7jVHKm9k1q3K/d1INbhN8sI2wqdh+hunmFje86dJaz0xrzRknmnV5xe4+It+cjJtIO9pTkMrr405KcTkt6DXdCfv1i+52L5E3vV73UbYWEHVGKVlzDcP/6qTbfOOdUZw8Tb82XxsouZ1dJS4fm2rfPMbJ2J5g+Cz2/DuWn8G7J9SGhIt+dlt9Jr5Vs8voJ7li6hJtNaJPpBI9oH5llPR9y9KYnJb0Gu6k3dr2yqWNtvbcH5TqFLVYqVlLNyPtTK9DaaWs2+Qljeu+PXjoslbl8hrv2FUUda80pB53L6fFhZIjmTHDXsbHjx+aZ/qfTqzjYlZMSa3Bb3Gg7xxIyzTz/v8trjCVnlLe1GrblMi7+PbMG1Sq/cp+ezoFkWjA7pYXlOg+SOMeaUhWSk3F7DP9D5zPuZTnTeqUfBZspuTonTBuZj8FvQaD/LGn2jzFquvp+dlK2yVN224OmyTXfXji8C0ec37ozYpI3a00CJo8M7F8pqwrLylIdm/cstW6LKadFjcf1AHYao0nZZE6oJzMfkt6DUe2rx5efWByLRkha1t3gZ9Pan3puThyKuIDi2mibx2R+gfkLdAtqdgYb6mNZvyJt27e3ob4jJakTduU2blzRbbmbwlIXl5S7fCXd5sgfa4dF/N2/sDtZgO5c22TWcZLwvyRp/1S/NCZreZJ2VltNFsKJE3r0ZH8paG5Nu8+QXs9rq3ebNd7Km8m23eTEwLZ1MOxWHOsL3PHrunB2zpqR+zv+a7zrAlZTQ4kMg1G8Z5eZMDcPt+diJveYj9WJZthTGvSW9DulOjSmK80dsQFWl6GzIxuS3oNV2ObYgE1P0IG2Mb0q6yjEnxYffm2IaSMhp1ldlI26WZk9d0faqW4bibmvcP4ZYQ3c+r+p6nhQVMV/FTvaAtoiaqn/fJz2ps2dPM6eFCP28ak9uCXtPpqLJlfjyTpXCSIu0SXb9ODhbiUWXFMmZNTlIk5cenllJ5kzNLT8IO5I02W3/8SkPm5jX9ocwv8OCh2Vq9x2ZBYXxeBVSVYDh/m25O/gxbLia7Bb2m2/G8ZiTp5nje4unhdEhOcTzvZhnra/eD4HQ87+P4pH5GXnNOXw9F6EDe9QtjZ2lIYWxDssAv30d/vB19HBdRFz0iOnj0Mne8e6d2a7IbsjHZLeg1PmZS2IE5jjRoNYjjQM3O/rduy/EyDWje5Hi2yZBIaSDvTrzI26TqbTQYXRrIuxM/EzAbzERrNA1IGsi7E0+zh50VbDYBUxrIuxOmvoNYkBfEgrwgFuQFsSBvkYU+CWcODM05ub6f5B8oyFvEyDsu/Aw9A3mLGGHj0YrI21uQt8giHWFlR2Qhby9B3iJK3t+ZQZrL0cnvY3l/+jId3LY2g9/MODA9COvEXGFmOTJDk1Vb2VwoZPz3kR6LmF0dWgF5iyh5/32i2w3z4PRbK+88Mw7YtibMAMVkCr+5ANumvL8NshfFORnSCcKuQd4iSt6xHmGsZiG9toPIk1Hs9lI3t8pHs8zJcz2M3l5pZENexWf51aEdkLeIEm2qvvD1rC+tsm78XuiZMmaiWXyB4alaxly/Ln5tU95LM4Mnszq0A/IW0fKqmVy3amaokTd6ambPZC5LbQTVPcHn/22nfhTl1auVrQ5NQd4iWt7IvtObWfTPyJt87dtLcax//K+Htp1r279mEk1R3ovbMCysDm2AvEW0vMrJbyam9s3Lq9q5L+In08wERzsrV1ex2+Wl0dsayFtkoQ2LvHswstrq/9MqU+v64D//NrEnkXVfmel7iOU17eGsvNS4rYO8RYy8pgdsqh/HsYwG+2Q1SV+7f2OuY5dp3Gblza0OLYG8RYy8+kDM3nHL9DacvtUPl3E9aoyMK9VFLG/wb/oCIjl5c6tDSyBvESOvPhCzOm7085qX9E22pkZIfVlxczGltHWbkZd+3i5A3iKL5EyE0tHKG59v0D/PEkHHmTNs2UMzdQmarLy51aEdkLdIrlWg61JtnL5yjxnBsFYt3Ec/2F5bM9Dhkb3CjL74zbNfrjfkza4O7YC8IBbkBbEgL4gFeUEsyAtiQV4QC/KCWJAXxIK8IBbkBbEgL4gFeUEs/wJhcU7X3WHl5QAAAABJRU5ErkJggg==" />

<!-- rnb-plot-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


# Assess performance in external validation

II. Look at performance as a function of prevalence

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHJldmFsZW5jZXMgPSBzZXEoMC4wMSwgMC45OSwgMC4wMSlcbmJyaW5hdGkuc3luLmZhY3RvcnkgPSBnZW5lcmF0ZVN5bnRoZXRpY0RhdGFGYWN0b3J5KGJhc2VERiA9IGJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCwgbWV0aG9kID0gJ2Jvb3QnKVxuZGZzX3ByZXZzIDwtIGxhcHBseShwcmV2YWxlbmNlcywgZnVuY3Rpb24gKHByZXYpIGJyaW5hdGkuc3luLmZhY3RvcnkocHJldiA9IHByZXYpKVxuZGZzLnByZXZzLnBlcmYgPC0gbGFwcGx5KGRmc19wcmV2cywgZnVuY3Rpb24oZGYpIHtcbiAgYXNzZXNzUGVyZihwcmVkaWN0ZWQgPSBwbG9naXMocHJlZGljdChnbG0ub3JpZy5maXQsIG5ld2RhdGE9YXMuZGF0YS5mcmFtZShkZikpKSwgb2JzZXJ2ZWQgPSBkZiRTV0FCKVxufSlcbmBgYFxuYGBgIn0= -->

```r
```r
prevalences = seq(0.01, 0.99, 0.01)
brinati.syn.factory = generateSyntheticDataFactory(baseDF = brinati_covid_study_v2.imputed, method = 'boot')
dfs_prevs <- lapply(prevalences, function (prev) brinati.syn.factory(prev = prev))
dfs.prevs.perf <- lapply(dfs_prevs, function(df) {
  assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(df))), observed = df$SWAB)
})

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Calculate the different validation measures

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuXG5kZnMucHJldnMucGVyZi5kZiA8LSBkYXRhLmZyYW1lKG1hdHJpeCh1bmxpc3QoZGZzLnByZXZzLnBlcmYpLCBucm93ID0gbGVuZ3RoKHByZXZhbGVuY2VzKSwgYnlyb3cgPSBUUlVFKSlcbmNvbG5hbWVzKGRmcy5wcmV2cy5wZXJmLmRmKSA8LSBuYW1lcyhkZnMucHJldnMucGVyZltbMV1dKVxuZGZzLnByZXZzLnBlcmYuZGYkUHJldmFsZW5jZSA9IHByZXZhbGVuY2VzXG5kZnMucHJldnMucGVyZi5kZiA8LSBnYXRoZXIoZGZzLnByZXZzLnBlcmYuZGYsIE1lYXN1cmUsIFZhbHVlLCBEeHk6YFM6cGAsIGZhY3Rvcl9rZXkgPSBUUlVFKVxuYGBgXG5gYGAifQ== -->

```r
```r

dfs.prevs.perf.df <- data.frame(matrix(unlist(dfs.prevs.perf), nrow = length(prevalences), byrow = TRUE))
colnames(dfs.prevs.perf.df) <- names(dfs.prevs.perf[[1]])
dfs.prevs.perf.df$Prevalence = prevalences
dfs.prevs.perf.df <- gather(dfs.prevs.perf.df, Measure, Value, Dxy:`S:p`, factor_key = TRUE)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Plot the distribution of validation measures

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxubWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpXG5wPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgK1xuICBnZW9tX3BvaW50KGRhdGE9ZGF0YS5mcmFtZShNZWFzdXJlID0gbWVhc3VyZXMsIFZhbHVlID0gYygxLCAwLCAwLCAxKSksIHNpemU9MywgY29sb3I9J2JsdWUnKSArIHRoZW1lX2J3KClcbnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApXG5TYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1ib3hwbG90LnBuZycpXG5wXG5wPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1QcmV2YWxlbmNlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSwgZ3JvdXA9TWVhc3VyZSkpICsgZ2VvbV9saW5lKHNpemU9MS4zKSArIHRoZW1lX2J3KClcbnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApXG5TYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1saW5lYXJwbG90LnBuZycpXG5gYGBcbmBgYCJ9 -->

```r
```r
measures = c('C (ROC)', 'Brier', 'Intercept', 'Slope')
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Measure, y=Value,  color=Measure)) + geom_boxplot() +
  geom_point(data=data.frame(Measure = measures, Value = c(1, 0, 0, 1)), size=3, color='blue') + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-boxplot.png')
p
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Prevalence, y=Value,  color=Measure, group=Measure)) + geom_line(size=1.3) + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-linearplot.png')

<!-- rnb-source-end -->

<!-- rnb-plot-begin -->

<img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArwAAAGwCAMAAAB8TkaXAAABAlBMVEUAAAAAADoAAGYAAP8AOjoAOmYAOpAAZpAAZrYAv8QzMzM6AAA6ADo6AGY6OgA6Ojo6OmY6OpA6ZmY6ZpA6ZrY6kJA6kLY6kNtmAABmADpmAGZmOgBmOjpmOmZmOpBmZjpmZmZmkJBmkLZmkNtmtrZmtttmtv98rgCQOgCQOjqQOmaQZgCQZjqQZmaQkDqQkGaQkLaQttuQtv+Q27aQ2/+2ZgC2Zjq2kDq2kGa2tra2ttu225C227a229u22/+2/9u2///HfP/bkDrbkGbbtmbbtpDbtrbb27bb29vb2//b/7bb/9vb///r6+v4dm3/tmb/25D/27b/29v//7b//9v////uJwYKAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAZKklEQVR4nO2dDXvbxpVGQdnSSmGTOLU3YuJsmnVjVt1qk3XXtLPdRpW5duI0akmZxP//K4v5wBcBUcDgkhckz3keWyIE4oXIo+HFYAaIYoAdJdLeAYBQkBd2FuSFnQV5YWdBXthZkBd2lm7yfgSQIaRkczrKK7QXCf+Q2xThKuHIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIq4Nm+NnZmWI68kpwqPKenanaezjyLl8No+j4Rb4AeTuDvOG0kXeeqGs4z5Ygb2eQN5wW8i5G0fFVvHwdDV6mi5C3O9S8wbSQdxqd3Jivk7zpRV7CM/os7/IiGq8uQ17CM/os72J0dL26DHkJz+izvPPhyc27p/Q2EH4HfZZ3Fp1cFnsb7GznfwB4+i1vFD25iZevorz2peUlPKPn8romdxKdpsuQl/CMfsvr+3fnw+zIDXkJz+izvJmzyEt4HX2WdzHyte7Mn6yIkZfwAn2WN6t1qXkJr6PX8s6H0eMr09vA2AbCa+i1vPHMjyrLzxIjL+EZ/ZY3vr1M9H10lS9AXsIzei5vBeQlPAN5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IFwB5dSBcAOTVgXABkFcHwgVAXh0IF2Dn5NW+ZSj0h52TV2gv4v1pfw43HHl1IFwA5NWBcAGQVwfCBUBeHTTDz87OFNORV4JDlffsTNVe5JUAeVVAXgmQVwXkleBQ5aXmDQd5CRcCeXUgXADk1YFwAZBXB8IFQF4dCBcAeXUgXADk1YFwAZBXB8IFQF4dCBcAeXUgXADk1YFwAZBXB9VwzaENyCsC8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8mrw5s0bsW0hrw6HKu+bN4L2Iq8OyCtAz+V99zSKBs9u8gXIKwDyhtJG3mlkOb7OliCvANS8obSQdz4cPI/j21F0mi1CXgHobQilhbyT6Nx8mQ+PsqYXeQVA3lDaH7AtRsgrCvKG0l7eWXSSHbIhrwDIG0pred8Oo3H2AHkFQN5QWso7iaLBn923Hxm073e7D5xp74AQfZd3efnpMDJ9Dh5aXgFoeUNpX/O+K9QNyCsA8oYScHp4mh+xIa8AyBtKgLyFjl7kFQB5Q2ku7/LClwvIKwvyhtLqDNtp6WuMvCIgbyitxjZET27i5ato8DJdhLwCIG8obWret25UGV1lsiBvKK0O2G6/StR9fJUvQF4BkDcUZlKohyNvKMirHo68oSCvejjyhoK86uHIGwryqocjbyjIqx6OvKEgr3o48oaCvOrh+y2vHwmzGJ3c1P68C8irHr7v8trRBPMh8m6Ig5V38xcdmQ+/MVdMmH6DvBviUOXdwuWe5sNvv76Jl3/8PpF3eRFFpoiYD6PIXATEfh3H898kbXPy3/zj75Ifpys1AHnVw/dd3uf/cR3PP3t7crO8OE2a4JObxWicfD26dtIOx5m8w1MzbNyt1CQRedXD913e8XQcz85nJzcz06Am5n4wZtqW1jWwubyJ1OlKTRKRVz1832ve8ew0nowTef2FGs/NhWsiMyp8Etl5DZm85mu+0v0gr3r4vvc2jBdf//3LayOvLwYWo0Rca2rybZSWD6m8LQ7skFc9fO/ljV99fxqbssFPwZkZP9MHSYVQkneWz9O5F+RVD99/eadJFTAzB2yJtYmcxk/T/WvL28TYxeg8OU4bOHnTlZokIq96+P7La6yc+a4yo6W5aNh/J8dkvvS1XWbffuHkzVZqAPKqh++3vJsEedXD97u3YZMgr3r4fvfzbhLkVQ9H3lCQVz0ceUNBXvVw5A0FedXDOWALBXnVw+kqCwV51cNpeUNBXvVwat5QqvJ++Okvjcf1IK8AByPvP6WSUlblfTs0kzAWX/y+0bORVwDkDWVF3ldmIHAib/Hu2GtAXgGQN5SyvLMokzdqMhEDeQXYd3n/WUUosSyvmZfxd3Nj7GkUNRnRjrwC7HtvQ1XVjcibNLiDl/au7nZ6xv3P5vatAuz77Vu3J68pGVrIK7QXMS2vEls4SbFeXlOpph/y03N7JQen3jTypatfY/5ZpRS4o+WdUTZsjYOWd2pEnTjXEj3dbHgzCdP8M/OD8jWmlRnFlZr35K+jo7+ZDrMm3Q3IK4CqvHIlb5C8iy/czB/bxE7GsZM3+X/xZfpNtoZbVKQsr2u13Q2rmkwjQl4B9l7eNb0Ns8Lnu5Eza3ndDxJnC2tMVjvAVvp532byNrrqA/IKsPfyrlkyK3y+m6LWtZ5m/rC/2+q4sEalblg9w2ZutZZw/KLT/gaAvCooy1s8CjNtrGl57bx352zS8hbWuFfehPfv33fd3wCQV4Ue1by25TVlgzlEq6l5G8krsL8BIK8KyvIWexuymtddKNJ22BZ7G+6teYX2NwDkVUFb3mI/b97bMPQX4yv2897X22CGNGRs9yTFgwcPxLbVnj2S900rkpe93RPWJHce21BzGqLIPf28avI+cAhtrT37JG+Lddu/7CHyZtx7VriqZ4EGZ9iQd+sgbygleZc//Y/jMoqeNHm2kLwPHijbu0/ytqoZHJstGzI2PpPCM93qeF7klaPH8opzh7xJAbHFgTl7Ju/ZZlkf3qJsCHjZd0Veat5QNtv3JSfvtmtecerlXU66H7DV9JDciX8V2zwl7Ne9g32St33d0OYpa5KV5S31NnQtG9qYaPVt94SOv3eZPZK3HZIfd/fK+wexKM/d8nY9YGvnYmu6/dorHKy82zjDlrE1eQeNLtwg9kkh72MrkFcAZXmzft6/NBxXhrx1IG+JP1QRSuzJwBzkVdr6NuRtsCQM5DUgrwAHK2+s6i7yShAir5v24/oG0vtoNyeTt9TTsP0hkfs0nneX5N3C1Pe18n587W+cHQLyyocjb4l75V3+8eX80dOj//04vdiNedSkFUZe+XDkLbGut8G1vJ/d2Ju8Jt9PzuPpqb3laxPymvfD+zK/dNjfAJBXZevaLa+f6m7E9bMtF19eN61++3LAhrw6W9eW1024HKfyjuzlbrrL+6HBs5G3DuQtca+88SST182xDJX3px++cTQqmZG3DuQtca+8SbHg5TU17+zkJkze5QUHbN1B3hJN+nlTee11SgNb3im9DQIgb4ktjW0wDe/gkyh6+GkUHTeRH3nrQN56Nj4kcuz+W064oUowyFvPpuW114UylzedD2tmUtxeJuXE46t8AfLWgbz1bFpec3XJaWSvqF6tef1JuMJlp/dFXtmXFXm3xGrNe/zCXNjs+fJ19YAt+enJVXx7UZjd1lt5aw4SJFkfjrxbotLb4O4gWDcBcz60Ptua2LMn8j548EBP3rZzIJE3o3L71kTQiZW3ctGzqb/HyjT/0V7I6yeA68jb/tIJyJuxeobt56+v46U5LqteqyztgJjljXJv5W2D+BVPkLce8XLqjouOvK92NSQlrx/xPtwreeWvNdXiTQoIR96MfDzvv71Yf5UR1xURZ7XvRwbtW4YKkPkjtsUNyyu2n8Joypscpj1adw+gVXkNtLx10PKWaH2twMaszqR4/ONda+6rvBuoeVvgw+WuEtmSbcjbYEkY1WlAg+IptAL7WvMirxCK8sbx7ffpzVsHz+rmAO1rb4P43Vzo5y2xTt588OP9wyCra5R7G3J/Hz6rHL7tUj+vajhn2EpsS96EDz8M7xjPu0tn2FTDkbfEPfLOP/tTFI3tnPfCxHc3Kt0tWHz5XfKg5nrndf28zt+6sQ3HuzK2QTUceUusK9uNvMNs6k8+8d3crjVZ6BYsRic3s6OaOcV18i5/+KR2JoWbs7GPo8qQV4DQltfNACpPfHei+gWL0dhemKRL2cB43oYgb4nm8uYT390dA/0Ce+9sN8O4zF0HbGZsZPj+BoC8Klvvlbz5xHff8roFRt77Wt7by4K5DW5IsWZ/A0Bela33Sd7CxHdT8yZf3YLF6LR+QnzNtcqamrtmfwNAXpWt90be5UXaxeCu5ZD2Npj/v/id+WLWKG9nVd6a/t2A/Q0AeVW2vgtjG2zNW0dJ3nbmxshbD/LWE7rnDeT9qq25MfLWg7z1bGkwelOQtw7krQd5NwLyCqA/h60dyFsH8m4J5JUPR94tgbzy4ci7JZBXPhx56xG855sDeeXDkbce5N0IyCsA8uqAvALcJcObKkKJyCsfjrwlqqoiryTIKwDy6oC8AgTJa8aDHV03vv9PGeSVD0feEuvkNdPT4mnzO6+VQV75cOQtsU7e4k0Ei/Pc47hmpnsF5JUPR94S63oblhfuOgr5JCA/z91Pg18P8sqHt7nwWACi+6pd85o7oLipP2a2ZdIG+3nuftb7+kTkVQ/fbEN9D+ryJsx/89LPZfPWxpOxn/W+PhF51cMPWd6ZrQzsNRnSltfNc7+31TUgr3r4Ictrexv8AVtpnrt/tD4RedXDD1le28/rp7v73gY3z93Pel8P8qqH7728rcY23DlVuAbkVQ/fd3kzGp0VRt6dCkfeUJBXPRx5Q0Fe9fCDkVcc5FUPR95QkFc9HHlD6Siv9i1D94He3o+1JTsnr9BexLS8Shxyyyu0FzHyKoG8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qqAvBIgrwrIKwHyqoC8EiCvCsgrAfKqgLwSIK8KyCsB8qpwOPK+expFg2eFW3EjrwDIG0obeaeR5fg6W4K8AiBvKC3knQ8Hz+P4dhSdZouQVwDkDaWFvJPo3HyZD4+yphd5BUDeUNofsC1GyCsK8obSXt5ZdJIdsiGvAMgbSmt53w6jcfYAeQVA3lBayjuJosGf3bcfGbTvd7sPcO/hUNrJu7z8dBiZPgcPLa8AtLyhtK953xXqBuQVAHlDuV/emTs1kRe60/yIDXkFQN5QAuQtdPQirwDIG0rzsmF54QVGXlmQN5RWZ9hOS19j5BUBeUNpNbYhenITL19Fg5fpIuQVAHlDadPb8NZVv3SVyYK8obTqKrv9KlH38VW+AHkFQN5QmEmhHo68oSCvejjyhoK86uHIGwryqocjbyjIqx6OvKEgr3o48oaCvOrhyBsK8qqHI28oyKsejryhIC/hQiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuADIqwPhAiCvDoQLgLw6EC4A8upAuAA7J6/2LUOhP+ycvEJ7Ee9P+3O44cirA+ECIK8OhAuAvDoQLgDy6kC4AMirA+ECIK8OhAuAvDoQLgDy6kC4AMirA+ECIK8OhAuAvDoQLgDy6kC4AMirA+ECIK8OhAuAvDoQLgDy6qAZfnameQtM5JXgUOU9O1O1F3klQF4VkFcC5FUBeSU4VHmpecNBXsKFQF4dCBcAeXUgXADk1YFwAZBXB8IFQF4dCBcAeXUgXADk1YFwAZBXB8IFQF4dCBcAeXUgXADk1YFwAZBXB8IF2Dl5ATKElGxON3kF2f6vTrh2eFeQl/CdBXkJ31l6Iy9AW5AXdhbkhZ0FeWFnQd6tsbzpwyb2iQ3L+/PlMIoGj69Wl0+OruN58iPHwyfZm/Lhh0/Mgme/3LGNWTTusDuLUZr4rKTBNDrtsNUK82Hy663y7mmXPe+yieW7r8yv7F/BxWjwsuuO9ISNynv7NNXzvPwDq2AubxT5N3v5KluS+ry6jUmNF43J5I2ik6K9W5B3edHpz67DJhbZK/jEPkTeJiSqHL9IHLm9XLF3MTLq5G/xu6GTJ3l3osc/xqaxTe2qbGM+7ODZYuQj08TN0Ct5k2cdm9f0w2v3CiJvE5JXLW3eplHp3ZzaN6HwFvtvJ1H2us6dXTXbmHR47TN5k7a/Swt+D72SN98Z9zsjbxMKfiQve6HpdQ1v8S12L+h8WHhZZ1bkmm10aXpzed13s+h8NoyOX/qy4dZU14+uXLz7SVCM+83mw5Nff0g2ePwiNn96hvFdIe+Sj/aHz92zSyu4TdhyNd9EK2bZn797lVN5zWdZXUx5D3rNBuWdFIT9UK4wvYWZmM7RcuVpn16zjeVFeMuRy5uYdWNyfzs0BbdLnvkifBznPwmKSeU9fppV65l5tSHTQlVfXuHkT/bBYNxB3qOShl5eHzOoxJT3oNdsTt47P+RS+zJ5l29tjbDyBNNi1G5junr415xizXtuM6KTq/gX92eTlNfJUeLytW/y3U+CSOWNos9v4tsL+zfgf5XakGTNz2+SZe5TvbyCeVTeREuSZw3sgUT2GrhPORP94bImprQHvWZz8t5ZW6XSFnsbjq8rTapZrXYbs/BjrUJvg2vl/Btk5U0bfvt11uWty+T1nzBmU9682pBpWt+PKytkR7LnHXsbHj5+4R/a13Tii4lJNaa0B71GQd60CCv08z6/qT7hTnlre1Gb7lMm7+ObOC+pzfuU/e3YiqLTAV0qr9ug+yWcebUhRSlXV/CP7GsWfMxnOm9MUfB59jJnm7IbLsWU96DXKMib/kW7t9h8PD2ve8Kd8uaFa8A++ae+v4xczeveH7NLBbGTlWZRh3culdeFFeWtDSn+lnfshd1Wlw6LD+/MQZjZmk3LIu2GSzHlPeg1CjVvWV57IDKuecKdNW+Hvp7ce7flw5HXkBxajDN5/Qthv0HeCsWegpn7mLasypt1797T25BuQ0TetKYsylvc7MbkrQkpy1u7F+HyFjfoj0vva3l7f6CWskF5i7XppOBlRd7kb/3ULSi8bO5B3TYkyoYaectqbEje2pByzVtewe9veM1b7GLP5V2teQsxAmdTtsV2zrC9LR675wds+akf93pN151hy7bR4UCiVDacl+XNDsD9+7kReetD/J9l3V4487r0NuQvatJInK/0NiSbdL0NhZjSHvSaTY5tSAS0/QgrYxvyrrKCSelh9+rYhpptdOoq85G+S7Mkr+v6NJXh+WZa3t/Hd4TYfl7T9zyurOC6ip/aFf0mWmL6eZ/8YsaWPY3y08OVft48prQHvWajo8rm5fFMnspJirxLdPk6O1hIR5VVtzHpcpIi2356aimXNzuz9CTegLzJbts/v9oQf4bN/lGWV3j4idvbl3G+iZaYJsFxfJXvTvkMWymmuAe9ZrPjed1I0tXxvNXTw/mQnOp43tVtLC/CD4Lz8byP05P6BXndOX07FGED8i4vnZ21IZWxDdkKv75Kfnk/+jjdRFvsiOjo0YvS8e6teVmzl6EYU9yDXqMxk8IPzAmkQ9Wwc2yp7Ox/dVuPyjSgaZfj2S5DIncN5F2Lirxdmt5Og9F3DeRdi84EzA4z0TpNA9o1kHctSrOHgxXsNgFz10DetTD1HXYW5IWdBXlhZ0Fe2FmQt8rMnoRzB4bunFzfT/IfKMhbxcl7XvkeegbyVnHCpqMVkbe3IG+VWT7Cyo/IQt5egrxVjLy/dYM058PB71J5f/4qH9y2dIPf3DgwOwhr4K4wMx+6ocmmVnYXCjn/v6Edi1h8OoiAvFWMvP86snXDNDr6zss7LYwD9tWEG6CYTeF3F2BblffTqHhRnMEhnSDcNMhbxch7bkcYm1lIr/0g8mwUu7/UzY3x0a0zeG6H0fsrjazIa/i8/HSQAXmrGNHG5gPfzvqyKtvi98TOlHETzdILDI/NOu76demyVXlP3QyewtNBBuStYuU1M7luzMxQJ2/y0M2eKVyW2glqe4KP/8tP/ajKa59W93ToCvJWsfIm9h1dT5J/Tt7sY99fimP5039+4utcX/+6STRVeU9u4rjydJAAeatYeY2T345c61uW19S5l+mDcWGCo5+Va5vYu+Wl6BUDeavMrGGJdw+HXlv7f95kWl0f/vtfR/4ksu0rc30PqbyuHi7KS4srDvJWcfK6HrCx/XqeyujwDxajfNmHH9x17ArFbVHe0tNBCOSt4uS1B2L+jluut+Hoyn45TdtRZ2TaqM5SeaN/sRcQKclbejoIgbxVnLz2QMzruNLP6xbZm2yNnZD2suLuYkp5dVuQl37eTYC8VWbZmQijo5c3Pd9gv59kgp4XzrAVD83MJWiK8paeDjIgb5VSVWDbUmucvXKPG8GwNBXuox99r60b6PDIX2HGXvzm2a8XK/IWnw4yIC/sLMgLOwvyws6CvLCzIC/sLMgLOwvyws6CvLCzIC/sLMgLOwvyws6CvLCz/D+3rBiqsWo1PgAAAABJRU5ErkJggg==" />

<!-- rnb-plot-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucFxuYGBgXG5gYGAifQ== -->

```r
```r
p

```

Calculate Decision Curve

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KYGBge3IgbWVzc2FnZT1GQUxTRX0NCmxpYnJhcnkocmVhZHIpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShybXMpDQpsaWJyYXJ5KHN5bnRocG9wKQ0KbGlicmFyeSh0aWR5cikNCmxpYnJhcnkocGxvdGx5KQ0KbGlicmFyeShzalBsb3QpDQpsaWJyYXJ5KHBlcmZvcm1hbmNlKQ0KIyBpbnN0YWxsIHRoZSAnbWknIHBhY2thZ2UgaWYgeW91IGFyZSBwbGFubmluZyBvbiByZWRvaW5nIGltcHV0YXRpb24gb24gdGhlIG9yaWdpbmFsIGRhdGFzZXQuIEZvciBub3cgdGhlIGltcHV0ZWQgZGF0YXNldCBoYXMgYWxyZWFkeSBiZWVuIHNhdmVkDQojIGxpYnJhcnkobWkpDQpzb3VyY2UoJy4vaGVscGVyIGZ1bmN0aW9ucy5SJykNCmBgYA0KDQpSZWFkaW5nIGluIG9yaWdpbmFsIGRhdGFzZXQsIHBlcmZvcm1pbmcgbXVsdGlwbGUgaW1wdXRhdGlvbiBhbmQgd3JpdGluZyBvdXQgaW1wdXRlZCBkYXRhc2V0LiBUaGlzIGNodW5rIGNvbW1lbnRlZCBvdXQgYmVjYXVzZSB3ZSd2ZSBhbHJlYWR5IGRvbmUgdGhpcyBhbmQgc2F2ZWQgdGhlIHJlc3VsdGluZyBmaWxlIGFzIG91ciBzdGFydGluZyBwb2ludC4NCmBgYHtyIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQojIA0KIyBicmluYXRpX2NvdmlkX3N0dWR5X3YyIDwtIHJlYWRfY3N2KCJkYXRhL2JyaW5hdGktY292aWRfc3R1ZHlfdjIuY3N2IiwNCiMgY29sX3R5cGVzID0gY29scyhHRU5ERVIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIk0iLA0KIyAiRiIpKSwgU1dBQiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMCIsDQojICIxIikpKSkNCiMgYnJpbmF0aV9jb3ZpZF9zdHVkeV92MiA8LSBicmluYXRpX2NvdmlkX3N0dWR5X3YyICU+JSBtdXRhdGVfaWYoaXMubnVtZXJpYywgLmZ1bnMgPSBmdW5jdGlvbih4KSB7eCswLjAwMX0pDQojIGJyYWluYXRpLm1pIDwtIG1pc3NpbmdfZGF0YS5mcmFtZShhcy5kYXRhLmZyYW1lKGJyaW5hdGlfY292aWRfc3R1ZHlfdjIpLCBmYXZvcl9wb3NpdGl2ZSA9IFRSVUUpDQojIA0KIyANCiMgYnJhaW5hdGkubWk8LWNoYW5nZShicmFpbmF0aS5taSwgeT1jKCdMeW1waG9jeXRlcycsICdCYXNvcGhpbHMnLCAnTW9ub2N5dGVzJywnRW9zaW5vcGhpbHMnLCAnQmFzb3BoaWxzJyksIHdoYXQ9J3R5cGUnLCB0bz0ncG9zaXRpdmUtY29udGludW91cycpDQojIHNob3coYnJhaW5hdGkubWkpDQojIGltYWdlKGJyYWluYXRpLm1pKQ0KIyANCiMgb3B0aW9ucyhtYy5jb3JlcyA9IDQpDQojIGltcHV0YXRpb25zIDwtIG1pKGJyYWluYXRpLm1pLCBuLml0ZXIgPSA5MCwgbi5jaGFpbnMgPSA0LCBtYXgubWludXRlcyA9IDIwKSANCiMgc2hvdyhpbXB1dGF0aW9ucykNCiMgcm91bmQobWlwcGx5KGltcHV0YXRpb25zLCBtZWFuLCB0by5tYXRyaXggPSBUUlVFKSwgMykNCiMgUmhhdHMoaW1wdXRhdGlvbnMpDQojIGRmcyA8LSBjb21wbGV0ZShpbXB1dGF0aW9ucywgbT0yKQ0KIyBicmluYXRpX2NvdmlkX3YyLmltcHV0ZWQgPC0gZGZzW1sxXV1bLDE6MTZdDQojIHdyaXRlLmNzdihicmluYXRpX2NvdmlkX3YyLmltcHV0ZWQsIGZpbGU9Jy4vZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2Jywgcm93Lm5hbWVzID0gRkFMU0UpDQpgYGANCg0KUmVhZCBpbiBkYXRhc2V0IGZyb20gQnJpbmF0aSBwYXBlci4gRGF0YXNldCBhbHJlYWR5IGJlZW4gaW1wdXRlZCBiYXNlZCBvbiBjb2RlIGFib3ZlLiANCmBgYHtyfQ0KDQpicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQgPC0gcmVhZF9jc3YoImRhdGEvYnJpbmF0aS1jb3ZpZF9zdHVkeV92Ml9pbXB1dGVkLmNzdiIsDQogIGNvbF90eXBlcyA9IGNvbHMoR0VOREVSID0gY29sX2ZhY3RvcihsZXZlbHMgPSBjKCJNIiwNCiAgIkYiKSksIFNXQUIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIjAiLA0KICAiMSIpKSkpDQpgYGANCg0KDQpgYGB7cn0NCg0KZ2xtLm9yaWcuZml0PC1nbG0oU1dBQiB+IC4sIGJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCwgZmFtaWx5PSdiaW5vbWlhbCcpDQpnbG0ub3JpZy5maXQNCndyaXRlLmNzdihzdW1tYXJ5KGdsbS5vcmlnLmZpdCkkY29lZiwgZmlsZT0nb3V0cHV0L09SLW9yaWdpbmFsLW1vZGVsLmNzdicsIHJvdy5uYW1lcyA9IEZBTFNFKQ0KcHJlZGljdGVkID0gcGxvZ2lzKHByZWRpY3QoZ2xtLm9yaWcuZml0KSkNCm9ic2VydmVkID0gYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkJFNXQUINCmFzc2Vzc1BlcmYocHJlZGljdGVkLCBvYnNlcnZlZCkNCnA8LXBsb3RfbW9kZWwoZ2xtLm9yaWcuZml0KQ0KcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCArIHlsaW0oLjEsIDIuNSkgKSANCiAgDQpwDQpnZ3NhdmUoZmlsZW5hbWUgPSAnZmlncy9PUi1vcmlnaW5hbC1tb2RlbC5wbmcnLCBwbG90PXAsIHdpZHRoPTYsIGhlaWdodD02LCB1bml0cz0naW4nLCBkcGk9MzAwKQ0KYGBgDQoNCiMgRHJhdyBST0MgY3VydmUNCmBgYHtyfQ0KY21zLmN1dG9mZnMgPC0gbGFwcGx5KHNlcSgwLjAxLDAuOTksMC4wMSksIGZ1bmN0aW9uKGN1dG9mZikgew0KICByZXQ9Y29uZnVzaW9uTWF0cml4KHByZWRpY3RlZCwgb2JzZXJ2ZWQsIGN1dG9mZj1jdXRvZmYpDQogIHJldCRjdXRvZmYgPSBjdXRvZmYNCiAgcmV0DQogIH0pDQpmZjwtZGF0YS5mcmFtZSh0KHNhcHBseShjbXMuY3V0b2ZmcywgZnVuY3Rpb24oY2YpIGMoIkN1dG9mZiI9Y2YkY3V0b2ZmLCAiU2Vuc2l0aXZ0eSI9Y2Ykc2VucywgIlNwZWNpZmljaXR5Ij1jZiRzcGVjKSkpKQ0KY29sbmFtZXMoZmYpIDwtIGMoJ0N1dG9mZicsICdTZW5zaXRpdml0eScsJ1NwZWNpZmljaXR5JykNCg0KDQpjdXRvZmZzX3RvX3Bsb3QgPC0gYygwLjA5LCAwLjMsIDAuNiwgMC45KQ0KcCA8LSBBcHBseUZpZ3VyZVRoZW1lKGZmICU+JSBtdXRhdGUoZnByPTEtU3BlY2lmaWNpdHkpICU+JSBmaWx0ZXIoQ3V0b2ZmICVpbiUgY3V0b2Zmc190b19wbG90KSAlPiUNCiAgZ2dwbG90KC4sIGFlcyh4PWZwciwgeT1TZW5zaXRpdml0eSkpICsgDQogICAgZ2VvbV9wb2ludChzaXplPTMsIGNvbG9yPSdkYXJrYmx1ZScpICsgeGxpbSgwLDEpICsgeWxpbSgwLDEpKw0KICAgIHhsYWIoJ0ZhbHNlIFBvc2l0aXZlIFJhdGVcbigxLVNwZWNpZmljaXR5KScpICsgeWxhYignVHJ1ZSBQb3NpdGl2ZSBSYXRlXG4oU2Vuc2l0aXZpdHkpJykNCiAgKQ0KZ2dzYXZlKGZpbGVuYW1lID0gJ2ZpZ3MvZXhhbXBsZS1zZW5zLXNwZWMtb24tcm9jLnBuZycsIHBsb3Q9cCwgd2lkdGg9NiwgaGVpZ2h0PTYsIHVuaXRzPSdpbicsIGRwaT0zMDApDQoNCmZvciAoY3V0b2ZmIGluIGxhcHBseSgxOmxlbmd0aChjdXRvZmZzX3RvX3Bsb3QpLCBmdW5jdGlvbihpKSBjdXRvZmZzX3RvX3Bsb3RbMTppXSkpIHsNCiAgcCA8LSBBcHBseUZpZ3VyZVRoZW1lKGZmICU+JSBtdXRhdGUoZnByPTEtU3BlY2lmaWNpdHkpICU+JSBmaWx0ZXIoQ3V0b2ZmICVpbiUgY3V0b2ZmKSAlPiUNCiAgICBnZ3Bsb3QoLiwgYWVzKHg9ZnByLCB5PVNlbnNpdGl2aXR5KSkgKyANCiAgICAgIGdlb21fcG9pbnQoc2l6ZT0zLCBjb2xvcj0nZGFya2JsdWUnKSArIHhsaW0oMCwxKSArIHlsaW0oMCwxKSsNCiAgICAgIHhsYWIoJ0ZhbHNlIFBvc2l0aXZlIFJhdGVcbigxLVNwZWNpZmljaXR5KScpICsgeWxhYignVHJ1ZSBQb3NpdGl2ZSBSYXRlXG4oU2Vuc2l0aXZpdHkpJykNCiAgKQ0KICBnZ3NhdmUoZmlsZW5hbWUgPSBwYXN0ZTAoJ2ZpZ3MvZXhhbXBsZS1zZW5zLXNwZWMtb24tcm9jLScsIGN1dG9mZltsZW5ndGgoY3V0b2ZmKV0sICAnLnBuZycpLCANCiAgICAgICAgIHBsb3Q9cCwgd2lkdGg9NiwgaGVpZ2h0PTYsIHVuaXRzPSdpbicsIGRwaT0zMDApDQp9DQpwIDwtIEFwcGx5RmlndXJlVGhlbWUoDQogIGdncGxvdChmZiwgYWVzKHg9MS1TcGVjaWZpY2l0eSwgeT1TZW5zaXRpdml0eSkpICsgZ2VvbV9saW5lKCkrIA0KICAgIGdlb21fcmliYm9uKGFlcyh5bWluID0gMCwgeW1heD1TZW5zaXRpdml0eSksIGNvbG9yPU5BLCBmaWxsPSdibHVlJyxhbHBoYT0wLjMpICsgDQogICAgeGxpbSgwLDEpICsgeWxpbSgwLDEpKw0KICAgICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKSANCikNCmdnc2F2ZShmaWxlbmFtZSA9ICdmaWdzL2V4YW1wbGUtcm9jLnBuZycsIA0KICAgICAgICAgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkNCnAgPC0gZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lDQogIGdncGxvdCguLCBhZXMoeD1mcHIsIHk9U2Vuc2l0aXZpdHkpKSArIGdlb21fcG9pbnQoYWVzKGZyYW1lPUN1dG9mZikpICsNCiAgICB4bGFiKCdGYWxzZSBQb3NpdGl2ZSBSYXRlICgxLVNwZWNpZmljaXR5KScpICsgeWxhYignVHJ1ZSBQb3NpdGl2ZSBSYXRlXG4oU2Vuc2l0aXZpdHkpJykNCiAgDQogIA0KZ2dwbG90bHkocCkNCmBgYA0KDQojIERyYXcgY2FsaWJyYXRpb24NCmBgYHtyfQ0KY2FsX2JyZWFrcyA9IHNlcSgwLCAxLCAuMSkNCmNhbF9kZWNpbGVzIDwtIGRhdGEuZnJhbWUoUHJlZGljdGVkPXByZWRpY3RlZCwgT2JzZXJ2ZWQgPSBhcy5pbnRlZ2VyKGFzLmNoYXJhY3RlcihvYnNlcnZlZCkpKSAlPiUNCiAgbXV0YXRlKEJpbnMgPSBjdXQoUHJlZGljdGVkLCBicmVha3M9Y2FsX2JyZWFrcywgaW5jbHVkZS5sb3dlc3QgPSBUUlVFKSkNCmNhbF9kZWNpbGVzX3N1bW1hcnkgPC0gY2FsX2RlY2lsZXMgJT4lIGdyb3VwX2J5KEJpbnMpICU+JSBzdW1tYXJpc2Uobj1uKCksIEF2ZXJhZ2VfT2JzZXJ2ZWQgPSBtZWFuKE9ic2VydmVkKSwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgQXZlcmFnZV9QcmVkaWN0ZWQ9bWVhbihQcmVkaWN0ZWQpKSANCg0KDQpmb3IgKGN1dG9mZiBpbiBsYXBwbHkoMTpucm93KGNhbF9kZWNpbGVzX3N1bW1hcnkpLCBmdW5jdGlvbihpKSBjYWxfZGVjaWxlc19zdW1tYXJ5JEJpbnNbMTppXSkpIHsNCiAgcDwtQXBwbHlGaWd1cmVUaGVtZUNhbEN1cnZlUG9pbnRzKA0KICAgIGdncGxvdChjYWxfZGVjaWxlc19zdW1tYXJ5ICU+JSBmaWx0ZXIoQmlucyAlaW4lIGN1dG9mZiksIGFlcyh4PUF2ZXJhZ2VfUHJlZGljdGVkLCB5PUF2ZXJhZ2VfT2JzZXJ2ZWQpKSApDQogIHA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApDQogIFNhdmVTdGRTcXVhcmVGaWd1cmUocCwgZmlsZW5hbWUgPSBwYXN0ZTAoJ2ZpZ3MvZXhhbXBsZS1jYWwtY3VydmUtcG9pbnRzLScsIGN1dG9mZltsZW5ndGgoY3V0b2ZmKV0sICAnLnBuZycpKQ0KfQ0KDQpwPC1BcHBseUZpZ3VyZVRoZW1lQ2FsQ3VydmVQb2ludHMoZ2dwbG90KGNhbF9kZWNpbGVzX3N1bW1hcnkgLCBhZXMoeD1BdmVyYWdlX1ByZWRpY3RlZCwgeT1BdmVyYWdlX09ic2VydmVkKSkgKQ0KcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNCg0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9jYWwtY3VydmUtZGVjaWxlcy1hbGwucG5nJykNCg0KcDwtcCArIGdlb21fc21vb3RoKG1ldGhvZD0nbG0nLCBzZT1GQUxTRSkNClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3MvY2FsLWN1cnZlLWRlY2lsZXMtYWxsLXdpdGgtYmVzdGZpdC5wbmcnKQ0KcA0KICAgICAgICAgICAgDQpgYGANCg0KQ2FsY3VsYXRlIGhvc21lciBsZW1lc2hvdyBzdGF0aXN0aWNzDQpgYGB7cn0NCg0KaGwub3JnIDwtIHBlcmZvcm1hbmNlX2hvc21lcihnbG0ub3JpZy5maXQpDQpobC5vcmcNCmBgYA0KDQpDcmVhdGUgc21vb3RoZWQgY2FsIGN1cnZlIGZyb20gVmFuIGhvb3JkZSBldCBhbA0KYGBge3J9DQpwPC1DcmVhdGVTbW9vdGhlZENhbEN1cnZlUGxvdChwcmVkaWN0ZWQsIG9ic2VydmVkKQ0KcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3Mvc21vb3RoZWQtY2FsLWN1cnZlLW9yaWctZGF0YS5wbmcnKQ0KcA0KYGBgDQoNCkJvb3Qgc3RyYXAgcGVyZm9ybWFuY2UgaW4gdGhlIG9yaWdpbmFsIGRhdGFzZXQNCmBgYHtyIHdhcm5pbmc9RkFMU0UsIGVjaG89RkFMU0V9DQpkZiA8LSBhcy5kYXRhLmZyYW1lKGJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCkNCm5Cb290ID0gMTAwDQpzZXQuc2VlZCgxMzQyMzQ5KQ0KdHJhaW5JZHhCb290IDwtIGxhcHBseSgxOm5Cb290LCBmdW5jdGlvbihpKSBzYW1wbGUoMTpucm93KGRmKSwgc2l6ZT1ucm93KGRmKSwgcmVwbGFjZT1UUlVFKSkNCg0KYm9vdFBlcmYub3JpZ2RhdGEgPC0gbGFwcGx5KHRyYWluSWR4Qm9vdCwgDQogICAgICAgICAgICAgICAgICAgZnVuY3Rpb24odHJhaW5JZHgpIA0KICAgICAgICAgICAgICAgICAgICAgYXNzZXNzU2luZ2xlVHJhaW5UZXN0KHRyYWluSWR4LCAoMTpucm93KGRmKSlbLXRyYWluSWR4XSwgZGYsIGdsbS5tb2RlbCA9IFNXQUJ+IC4sIG91dGNvbWUudmFyLm5hbWUgPSAnU1dBQicpDQogICAgICAgICAgICAgICAgICAgKSANCg0KbWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpDQpicmluYXRpLm9yaWcuYm9vdC5wZXJmIDwtIGRhdGEuZnJhbWUobWF0cml4KA0KICB1bmxpc3QobGFwcGx5KGJvb3RQZXJmLm9yaWdkYXRhLCBmdW5jdGlvbihib290KSBib290JHRlc3QucGVyZikpLCANCiAgbnJvdyA9IG5Cb290LCBieXJvdyA9IFRSVUUpKQ0KY29sbmFtZXMoYnJpbmF0aS5vcmlnLmJvb3QucGVyZikgPC0gbmFtZXMoYm9vdFBlcmYub3JpZ2RhdGFbWzFdXSR0ZXN0LnBlcmYpDQoNCmJyaW5hdGkub3JpZy5ib290LnBlcmYgPC0gZ2F0aGVyKGJyaW5hdGkub3JpZy5ib290LnBlcmYsIE1lYXN1cmUsIFZhbHVlLCBEeHk6YFM6cGAsIGZhY3Rvcl9rZXkgPSBUUlVFKQ0KcDwtZ2dwbG90KGJyaW5hdGkub3JpZy5ib290LnBlcmYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgKw0KICBnZW9tX3BvaW50KGRhdGE9ZGF0YS5mcmFtZShNZWFzdXJlID0gbWVhc3VyZXMsIFZhbHVlID0gYygxLCAwLCAwLCAxKSksIHNpemU9MywgY29sb3I9J2JsdWUnKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWUocCkNClNhdmVTdGRTcXVhcmVGaWd1cmUocCwgJ2ZpZ3Mvb3JpZy1tb2RlbC1ib290c3RyYXAtZXZhbHVhdGlvbi5wbmcnKQ0KcA0KYGBgDQoNCiMgQXNzZXNzIHBlcmZvcm1hbmNlIGluIGV4dGVybmFsIHZhbGlkYXRpb24NCg0KSUkuIExvb2sgYXQgcGVyZm9ybWFuY2UgYXMgYSBmdW5jdGlvbiBvZiBwcmV2YWxlbmNlDQpgYGB7cn0NCnByZXZhbGVuY2VzID0gc2VxKDAuMDEsIDAuOTksIDAuMDEpDQpicmluYXRpLnN5bi5mYWN0b3J5ID0gZ2VuZXJhdGVTeW50aGV0aWNEYXRhRmFjdG9yeShiYXNlREYgPSBicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQsIG1ldGhvZCA9ICdib290JykNCmRmc19wcmV2cyA8LSBsYXBwbHkocHJldmFsZW5jZXMsIGZ1bmN0aW9uIChwcmV2KSBicmluYXRpLnN5bi5mYWN0b3J5KHByZXYgPSBwcmV2KSkNCmRmcy5wcmV2cy5wZXJmIDwtIGxhcHBseShkZnNfcHJldnMsIGZ1bmN0aW9uKGRmKSB7DQogIGFzc2Vzc1BlcmYocHJlZGljdGVkID0gcGxvZ2lzKHByZWRpY3QoZ2xtLm9yaWcuZml0LCBuZXdkYXRhPWFzLmRhdGEuZnJhbWUoZGYpKSksIG9ic2VydmVkID0gZGYkU1dBQikNCn0pDQpgYGANCg0KQ2FsY3VsYXRlIHRoZSBkaWZmZXJlbnQgdmFsaWRhdGlvbiBtZWFzdXJlcw0KYGBge3J9DQoNCmRmcy5wcmV2cy5wZXJmLmRmIDwtIGRhdGEuZnJhbWUobWF0cml4KHVubGlzdChkZnMucHJldnMucGVyZiksIG5yb3cgPSBsZW5ndGgocHJldmFsZW5jZXMpLCBieXJvdyA9IFRSVUUpKQ0KY29sbmFtZXMoZGZzLnByZXZzLnBlcmYuZGYpIDwtIG5hbWVzKGRmcy5wcmV2cy5wZXJmW1sxXV0pDQpkZnMucHJldnMucGVyZi5kZiRQcmV2YWxlbmNlID0gcHJldmFsZW5jZXMNCmRmcy5wcmV2cy5wZXJmLmRmIDwtIGdhdGhlcihkZnMucHJldnMucGVyZi5kZiwgTWVhc3VyZSwgVmFsdWUsIER4eTpgUzpwYCwgZmFjdG9yX2tleSA9IFRSVUUpDQpgYGANCg0KUGxvdCB0aGUgZGlzdHJpYnV0aW9uIG9mIHZhbGlkYXRpb24gbWVhc3VyZXMNCmBgYHtyfQ0KbWVhc3VyZXMgPSBjKCdDIChST0MpJywgJ0JyaWVyJywgJ0ludGVyY2VwdCcsICdTbG9wZScpDQpwPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1NZWFzdXJlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSkpICsgZ2VvbV9ib3hwbG90KCkgKw0KICBnZW9tX3BvaW50KGRhdGE9ZGF0YS5mcmFtZShNZWFzdXJlID0gbWVhc3VyZXMsIFZhbHVlID0gYygxLCAwLCAwLCAxKSksIHNpemU9MywgY29sb3I9J2JsdWUnKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApDQpTYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1ib3hwbG90LnBuZycpDQpwDQpwPC1nZ3Bsb3QoZGZzLnByZXZzLnBlcmYuZGYgJT4lIGZpbHRlcihNZWFzdXJlICVpbiUgbWVhc3VyZXMpLCBhZXMoeD1QcmV2YWxlbmNlLCB5PVZhbHVlLCAgY29sb3I9TWVhc3VyZSwgZ3JvdXA9TWVhc3VyZSkpICsgZ2VvbV9saW5lKHNpemU9MS4zKSArIHRoZW1lX2J3KCkNCnA8LUFwcGx5RmlndXJlVGhlbWVMYXJnZUZvbnRPbmx5KHApDQpTYXZlTWVkUmVjdEZpZ3VyZShwLCAnZmlncy9jaGFuZ2UtaW4tbWV0cmljcy1kdWUtdG8tcHJldmFsZW5jZS1saW5lYXJwbG90LnBuZycpDQpwDQpgYGANCg0KIyBDYWxjdWxhdGUgRGVjaXNpb24gQ3VydmUNCmBgYHtyfQ0KY21zLmN1dG9mZnMgPC0gbGFwcGx5KHNlcSgwLjAxLDAuOTksMC4wMSksIGZ1bmN0aW9uKGN1dG9mZikgew0KICByZXQ9Y29uZnVzaW9uTWF0cml4KHByZWRpY3RlZCwgb2JzZXJ2ZWQsIGN1dG9mZj1jdXRvZmYpDQogIHJldCRjdXRvZmYgPSBjdXRvZmYNCiAgcmV0DQogIH0pDQoNCmNhbGN1bGF0ZV9uZXRfYmVuZWZpdF9hdF90aHJlc2hvbGRfZm9yX3NlbnNpdGl2aXR5X2FuZF9zcGVjaWZpY2l0eSA8LSBmdW5jdGlvbihzZW5zLCBzcGVjLCB0aHJlc2hvbGQpIHsNCiAgc2VucyAtICgxLXNwZWMpKih0aHJlc2hvbGQvKDEtdGhyZXNob2xkKSkNCn0NCm5ldGJlbmVmaXRfYWNyb3NzX3RocmVzaG9sZHMgPSBzYXBwbHkoY21zLmN1dG9mZnMsIGZ1bmN0aW9uKGN1dG9mZl92YWx1ZXMpIGNhbGN1bGF0ZV9uZXRfYmVuZWZpdF9hdF90aHJlc2hvbGRfZm9yX3NlbnNpdGl2aXR5X2FuZF9zcGVjaWZpY2l0eShjdXRvZmZfdmFsdWVzJHNlbnMsIGN1dG9mZl92YWx1ZXMkc3BlYywgY3V0b2ZmX3ZhbHVlcyRjdXRvZmYpKQ0KbmV0YmVuZWZpdF92c190aHJlc2hvbGRzID0gZGF0YS5mcmFtZSh0aHJlc2hvbGRzID0gc2FwcGx5KGNtcy5jdXRvZmZzLCBmdW5jdGlvbihjdXRvZmZfdmFsdWVzKSBjdXRvZmZfdmFsdWVzJGN1dG9mZiksIG5ldGJlbmVmaXQgPSBuZXRiZW5lZml0X2Fjcm9zc190aHJlc2hvbGRzKQ0KZ2dwbG90KG5ldGJlbmVmaXRfdnNfdGhyZXNob2xkcywgYWVzKHg9dGhyZXNob2xkcywgeT1uZXRiZW5lZml0KSkgKyBnZW9tX2xpbmUoKSArIHlsYWIoIk5ldCBCZW5lZml0IikgKyB4bGFiKCJUaHJlc2hvbGQgdG8gQWN0IikgKyB0aGVtZV9idygpDQpgYGANCg0K